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ABSTRACT 

Realistic  analytical  description  of  the  local  behavior  and  of  the 
overall  performance  of  crossed-field  accelerators  is  obtained  by  com¬ 
puting  the  coupled  two-dimensional  distributions  of  current  density, 
plasma  properties,  and  fluid  velocity,  temperature  and  pressure,  over 
the  whole  length  of  the  channel.  The  analysis  includes  consideration  of 
electron  nonequilibrium,  thermal  and  concentration  diffusion,  com¬ 
pressible  turbulent  boundary  layers,  finite  reaction  rates,  and  electron 
energy  relaxation,  and  the  effect  of  each  of  these  mechanisms  is  investi¬ 
gated.  The  development  of  compressible,  turbulent,  magnetohydro  dy¬ 
namic  boundary  layers  on  the  electrode  walls  is  obtained  through  a  novel 
formulation, that  includes  a  transport  equation  for  the  Reynold's  stress, 
and  a  fast  numerical  method  of  solution.  Analytical  results  obtained  in 
this  study  gave  consistently  good  agreement  with  experimental  measure¬ 
ments.  Applications  to  seeded-nitrogen  and  seeded-air  accelerators 
indicate  limited  current  fringing  in  the  entrance  and  exit  regions  and 
rapid  approach  to  periodical  current  density,  electric  field,  and  plasma 
property  distributions  in  the  main  part  of  the  channel.  Solution  of  the 
complete  coupled  problem  for  the  Hirho  experiment  gave  correct  pre¬ 
diction  and  interpretation  of  the  Hall  field  and  showed  the  substantial 
effect  of  local  axial  current  density  even  when  there  is  no  net  axial 
current  leakage;  it  also  showed  that  the  electrical  performance  could  be 
improved  substantially  through  optimization  of  the  electrode  configura¬ 
tion.  This  work  has  demonstrated  the  importance  of  realistic  compu¬ 
tations  in  the  design  of  efficient  magnetohydrodynamic  channels. 
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I.  INTRODUCTION 

This  report  presents  the  results  of  the  third  phase  of  a  study  aimed 
at  providing  realistic  performance  analysis  of  crossed-field  accelerators. 

The  basic  thesis  that  gave  motivation  to,  and  was  a  posteriori 
justified  by,  this  work  is  that  realistic  description  of  the  overall  perform¬ 
ance  characteristics  of  magnetohydrodynamic  channels  can  be  obtained 
only  by  accurate  knowledge  of  the  local  behavior  at  every  point  in  the 
channel  and  rigorous  analytical  account  of  the  physical  mechanisms  that 
influence  the  flowing  plasma. 

The  first  two  phases  of  this  work  were  devoted  to  obtaining 
rigorous  solution  of  the  electrical  problem  (i.  e. ,  the  current  density, 
electric  field,  and  plasma  property  distributions)  in  the  main,  periodic 
part  of  a  multielectrode  channel.  This  was  successfully  accomplished 
(Refs,  i  and  2):  A  realistic  two-dimensional  model  for  the  current 
density  field  and  for  the  distributions  of  electron  temperature  and  plasma 
composition  was  formulated,  and  a  numerical  method  and  a  computer 
program  were  developed  to  compute  these  interdependent  field  distribu¬ 
tions.  The  results  have  been  reported  in  Refs.  1-5  and  include  analysis 
of  the  effects  of  thermal  and  concentration  diffusion,  electrode -wall 
boundary  layers,  finite  reaction  rates,  and  electron  energy  convection. 

The  third  phase,  the  results  of  which  are  described  in  this  report, 
had  the  twofold  objective  to  (1)  extend  the  two-dimensional  solution  of  the 
electrical  problem  to  cover  the  entire  length  of  a  magnetohydrodynamic 
channel,  including  the  entrance  and  exit  regions,  and  (2)  to  obtain  solu¬ 
tions  of  the  gasdynamic  problem  (i.  e.  ,  the  fluid  velocity,  temperature, 
and  pressure  distributions),  including  velocity  and  temperature  profile 
developments  in  the  electrode -wall  boundary  layers,  over  any  prescribed 
accelerator  length. 

Fig.  1  is  a  schematic  diagram  of  the  analytical  problem  that  has 
been  solved  in  this  study.  Two  independent  variables  are  considered:  x  in 
the  direction  of  the  overall  plasma  flow  and  y  in  the  direction  of  the 
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Fig.  1,  Schematic  diagram  of  the  analytical  problem  solved  in  this  study. 
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applied  electric  field.  The  magnetic  induction  B  is  assumed  to  be 
directed  in  the  positive  z-direction  and  to  be  a  known  function  of  x.-.  ... 
Three  classes  of  computed  fields  are  distinguished.  First,  the  electrical 
unknowns,  namely  the  current  density  field  J  and  the  electric  field  E 
at  every  point  in  the  gas.  Second,  the  unknowns  that  characterize  the 
state  of  the  plasma,  namely  the  number  density  nQ  of  each  component  ' 
a  ‘that  is  needed  for  this  purpose,  and  its  characteristic  temperature 
Tq.  (Plasma  transport  properties  are  then  calculable  by  the  methods  of 
Ref.  6. )  Finally,  the  gasdynamic  fields,  namely  the  gas  density  p  or 
pressure  p,  gas  velocity  U  and  static  gas  temperature  T.  If  values 
for  the  gasdynamic  distributions  are  assumed  and  the  first  two  classes 
of  unknowns  (which  are  coupled),  are  solved  for,  it  is  called  solution  of 
the  (still  nonlinear)  "electrical  part"  of  the  problem  (see  Fig.  1).  This  is 
what  was  done  in  Phases  I  and  II.  If  values  for  the  current  density  or  the 
electric  field  distribution  as  assumed  and  the  last  two  classes  of  unknowns 
are  solved  for,  this  is  called  solution  of  the  "gasdynamic  part"  of  the  prob¬ 
lem  (see  Fig.  1).  In  the  work  reported  herein  (Phase  III)  all  three  classes 
of  unknowns  are  solved  for,  over  any  prescribed  channel  length,  by  itera¬ 
tion  between  the  electrical  and  gasdynamic  parts  of  the  problem. 

i 

To  obtain  the  solution  to  the  coupled  "electrical  +  gasdynamic" 
problem  over  the  whole  channel,  five  subproblems  had  to  be  considered: 

(i)  solution  of  electrical  part  of  the  problem  in  the  entrance  region  of  the 
channel,  (2)  solution  of  the  electrical  part  of  the  problem  in  the  periodic 
part  of  the  channel,  (3)  solution  of  the  electrical  part  of  the  problem  in  the 
exit  region,  (4)  solution  of  the  gasdynamic  part  of  the  problem  in  the  core 
of  the  flow,  and  (5)  solution  of  the  gasdynamic  part  of  the  problem  in  the 
turbulent  boundary  layers  over  the  electrode  walls.  Subproblem  (2)  had 
been  solved  in  Phases  I  and  II  of  the  project  (see  Refs,  i  and  2).  Thus, 
this  report  will  describe  the  solution  of  the  other  four  subproblems. 

Their  formulation  will  be  presented  in  Section  II,  the  methods  of  solution 
in  Section  III,  while  Section  IV  will  present  results  for  various  test  cases. 
The  discussion  will  be  more  detailed  in  the  case  of  the  more  complex  of 
the  four  subproblems,  namely  the  computation  of  the  compressible  turbu- 
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lent  magnetohydrodynamic  boundary-layer  development  over  the  electrode 
walls,  for  which  a  novel  formulation  (Ref.  7)  has  been  used. 

: With  the  completion  of  Phase  III,  this  work  can  make  the  claim 
that  it  can  "model"  the  electrical  and  gasdynamic  performance  of  magneto¬ 
hydrodynamic  channels  realistically.  Applications  have  shown  that,  unlike 
previous, idealized  or  oversimplified  theories,  it  can  predict  experimental 
results  with  accuracy,  and  can  provide  the  experimentalist  and  develop¬ 
ment  engineer  with  truly  rigorous  and  practical  assistance.  Towards  this 
aim,  every  effort  has  been  made  to  maintain  maximum  flexibility,  and 
the  .analysis  can  treat  (1)  different  channel,  electrode,  and  insulator 
geometries,  including  "staggered"  electrodes,  (2)  different  operating  . 
modes  and  magnetic  field  distributions,  (3)  different  current  leakage 
..rates,  (4)  different  combinations  of  boundary  layer  and  core  flow  initial 
conditions,  (5)  different  operating  fluids  (monatomic,  diatomic,  or  mix¬ 
tures  of  gases),  .(6)  different  seed  materials  and  seed  ratios,  (7)  different 
electrode  and  wall  temperatures  and  cooling  rates,  (8)  different  wall 
.ablation  rates,  and  (9)  different  load  conditions. 
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II.  FORMULATION 

1.  THE  ELECTRICAL  PROBLEM 

The  formulation  of  the  electrical  problem  has  been  described  in 
detail  in  Refs.  1  and  2.  Unknowns  of  the  problem  are  (1)  the  current 
streamf unction  ^  (2)  the  electron  temperature  Tg,  and  (3)  the  number 
densities  nQ  (a  =  1  ...  N)  of  all  the  components  that  are  needed  to 
characterize  the  state  of  the  plasma.  The  equations  of  the  problem 
.(equal  in  number  to  the  above  unknowns)  are  (i)  the  streamfunction 
equation,  (2)  the  electron  energy  equation,  and  (3)  the  continuity  equations 
for  each  plasma  component  (Refs.  1  and  2).  These  equations  are  uniformly 
valid  over  the  whole  length  of  a  multielectrode  channel. 

To  solve  the  electrical  problem  over  the  whole  channel,  it  is  con¬ 
venient  to  divide  the  channel  into  three  separate  regions: 

(1)  The  entrance  region 

(2)  The  main  (periodic)  part 

(3)  The  exit  region 

In  the  main  part  of  the  channel,  the  geometrical  periodicity  im¬ 
posed  by  electrode  segmentation,  together  with  the  fact  that  gasdynamic 
properties  change  very  little  over  a  streamwise  distance  of  the  order  of 
one  electrode  period  L  (Ref.  1),  leads  to  a  "quasiperiodic"  distribution 
of  all  the  unknowns  of  the  electrical  problem.  The  term  nquasi-periodicw 
is  used  to  indicate  that,  although  periodicity  is  a  sufficiently  good 
approximation  to  the  actual  situation  to  allow  the  use  of  periodic  boundary 
conditions  for  each  electrode  segment  in  the  mathematical  formulation, 
it  does  not  describe  a  physical  fact,  and  that  actually  the  plasma  proper¬ 
ties  do  change  slowly  from  one  electrode  to  the  other,  as  do  the  gas- 
dynamic  variables  and  possibly  the  external  loading  conditions  and  magnetic 
induction. 

The  "quasiperiodic"  condition  is  obviously  approached  asymptoti¬ 
cally  away  from  the  entrance  and  exit  regions  of  the  channel.  In  the  case 
of  the  first  and  last  electrode  pairs,  there  is  not  even  geometrical 
periodicity,  and  this  geometric  disturbance  will  of  course  affect  the 
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solution  of  the  problem  for  several  electrode  pairs  downstream  and  up¬ 
stream,  respectively.  In  addition,  even  if  geometrical  periodicity  exists, 
it  takes  a  certain  distance  before  periodicity  of  the  plasma  properties 
can  be  realized,  owing  to  the  finite  rates  of  ionization  and  of  electron 
energy  relaxation. 

1.  1.  Main  (Periodic)  Part  of  Channel 

As  mentioned  above,  this  part  of  the  channel  is  called  periodic  only 
because  periodic  boundary  conditions  can  be  used  for  the  purpose  of  ob¬ 
taining  the  solution.  This  part  covers  most  of  the  length  of  a  powered 
magnetohydrodynamic  channel.  The  solution  of  the  electrical  problem  in 
this  part  of  the  device  has  been  obtained  in  Phases  I  and  II. 

1.  2.  Entrance  and  Exit  Regions 

The  treatment  of  the  boundary  conditions  for  the  entrance  and 
exit  regions  does  not  differ,  since  the  situation  for  these  two  regions  is 
analogous  (they  are  mirror  images  of  each  other),  and  therefore  the 
formulation  discussed  in  this  section  will  be  valid  for  both. 

For  these  two  regions,  periodic  conditions  cannot  be  employed. 

In  addition,  owing  to  the  elliptic  nature  of  the  streamfunction  equation, 
the  whole  entrance  (and  respectively  the  whole  exit)  region  will  have  to  be 
treated  as  the  geometry  of  the  problem,  and  they  may  have  to  include 
several  electrode  pairs  in  order  to  account  for  the  transition  region, 
i.  e.  ,  before  periodicity  is  realized.  Experience  with  the  solution  has 
shown  that  the  axial  extent  of  these  transition  regions  does  not  exceed 
5-8  electrode  pairs. 

The  geometry  appropriate  to  the  entrance  problem,  shown  in 
Fig.  2,  starts  well  upstream  of  the  first  powered  electrode  pair  and  extends 
downstream  until  the  condition  of  quasiperiodicity  is  satisfied.  Thus  the 
entrance  region  actually  consists  of  two  parts,  a  powered  section  and  an 
unpowered  one.  The  unpowered  section  is  bounded  in  the  y- direction  by 
insulating  walls  of  length  Lq  and/or  Kq  unpowered  electrode  pairs. 

The  powered  section  is  bounded  by  K-Kq  electrode  pairs.  The  con- 
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ductor  and  insulator  segments  of  each  electrode  pair  k  (k  =  1,  2, . .  . ,  K) 
have  lengths  f  ,  and  i  .  ,  respectively,  so  that  the  total  length  of 
each  electrode  pair  is  =  ^ c  k  +  ^ i  k'  ^  staining  the  solution,  the 
number  K  of  electrode  pairs  is  varied  until  the  actual  extent  of  the 
entrance  region  (as  defined  above)  is  ascertained. 

The  boundary  conditions  on  the  y -walls  reflect  the  alternation  of 
successive  conductor  and  insulator  segments.  The  form  of  these  con¬ 
ditions  is  identical  to  those  used  in  Phases  I  and  II  (Refs.  1  and  2),  namely 
St  =  const  along  the  surface  of  insulator  segments  and  Ex  =  0  along  the 
surface  of  conductor  segments.  The  total  current  1^  k  passed  in  the  posi¬ 
tive  y-direction  per  unit  width  of  each  conductor  k  will  be  equal  to  the 
difference  between  the  value  of  the  streamfunction  St  on  the  insulator 
segments  k  and  k-1.  Thus,  in  the  case  of  an  unpowered  electrode  k, 
the  value  of  the  streamfunction  on  the  k  and  k-1  insulator  segments 
will  be  the  same. 


The  boundary  conditions  at  the  upstream  plane  x^  =  0  of  the 
entrance  region  will  be  of  the  Dirichlet  type,  i.  e.  ,  the  values  of  the 
streamfunction  for  all  y  will  be  known.  They  will  be  defined  by  the 
assumed  distribution  of  net  axial  current  leakage  Ix  through  the  channel. 
For  example,  vanishing  leakage  would  correspond  to  constant  St(0,  y), 
while  uniform  current  leakage  density  would  be  represented  by  linear 


variation  of  StfO,  y)  from  a  value  SP  on  the  anode  wall  to  a  value  S& 

Si  c 

on  the  cathode  wall  (St  -  ^  =  I  ).  At  the  downstream  (periodic)  end 

E3-  C  X 

L,  of  the  entrance  region,  one  of  two  kinds  of  boundary  con- 

Q  O 

ditions  can  be  used:  ( i)  Dirichlet  type  conditions  obtained  by  previous 


solution  of  the  one -electrode  problem  at  that  station  or  (2)  the  condition 
of  periodic  fields  at  the  stations  x^  and  x^-  Both  kinds  of  conditions 

have  alternatively  been  used  for  solutions  obtained  in  this  study  with  very 
little  difference  in  the  overall  results. 


Finally,  note  that  this  formulation  can  treat  the  case  of  an 

MHD  generator  (rather  than  a  crossed-field  accelerator)  simply  by  using 

negative  values  for  I  ,  and  I  ,  i.  e. ,  by  reversing  the  position  of  the 

y,  K  x 

anode  and  cathode  walls. 
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2.  THE  GASDYNAMIC  PROBLEM 

The  gasdynamic  problem  is  solved  in  two  parts.  In  the  core  of  the 
flow  a  quasi-one -dimensional  solution  is  used  (a  two-dimensional  inviscid 
solution  for  the  core  has  also  been  coded),  while  a  two-dimensional 
numerical  solution  is  carried  out  for  the  compressible  magnetohydrody¬ 
namic  boundary  layer  development  on  the  electrode  walls. 


2.  1,  Quasi-one -dimensional  Solution  for  Core  Flow 

This  formulation  assumes  that  all  quantities  represent  an  average 
over  the  cross  sectional  area  A  of  the  core  at  each  point  x,  so  that  only 
their  x-variation  is  of  importance.  For  example,  the  gas  velocity  u(x)  is 
defined  by  u(x)  =  u(x,  y)  dA. 

The  term  quasi-one -dimensional  is  used  to  indicate  that,  although 
one -dimensional  equations  cannot  by  definition  describe  fluxes  in  the  y- 
direction,  such  as  viscous  friction  and  heat  conduction,  nevertheless 
empirical  terms  describing  these  fluxes  have  been  included  in  the  momen¬ 
tum  and  energy  equations. 


It  is  assumed  that  the  area  variation  A(x)  of  the  channel,  the 
magnetic  induction  variation  B(x),  and  the  loading  variation  J^(x)  are 
given  over  the  entire  channel.  Note  that  >  0  in  an  accelerator. 


Under  these  assumptions,  (see  also  Ref.  8)  the  following  "channel 
flow"  equations  are  obtained  for  the  N  +  3  unknowns,  which  are  the  N 
number  densities  nQ  of  the  plasma  components,  the  axial  velocity  u, 
the  gas  temperature  T  and  the  electron  temperature  T  . 


(1)  Continuity  equations  for  each  component  a 

pu  4?  (na/p)  =  Sa  (1) 

where  sq  is  the  source  term  (i,  e.  ,  the  rate  of  production  of  particles  of 
component  a)  owing  to  the  reactions  that  take  place  in  the  plasma. 
Appropriate  expressions  for  these  source  terms  are  given  in  Ref.  2;  for 
example,  the  important  3 -body  recombination-ionization  reaction  of  a 
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species  contributes  to  the  rate  of  production  of  electrons  the  amount 


,  ,  ,  2 
a  =k*  nn-knn. 
r  r  e  a  r  e  1 


<!') 


where  k*  and  k  are  the  ionization  and  recombination  rate  coefficients 
r  r 

for  this  reaction  r  (they  are  functions  of  the  electron  temperature  —  see 
Ref.  2)  and  n.,  n  are  the  number  densities  of  ions  and  neutrals,  respec- 
tively,  of  the  species  in  question. 


(2)  Overall  momentum  equation 


du  dp 
pu  =  “ 


J  B  -  s 

y 


M 


(2) 


where  s^.  is  loss  of  momentum  due  to  friction  and  is  given  by 

SM  =  cfPu2/R 


(2‘) 


where  c^  is  the  friction  coefficient  and  R  is  the  hydraulic  radius  of  the 
channel.  In  the  case  of  a  channel  with  rectangular  cross  section,  R  is 
defined  as  R  =  aD/(a+D),  where  D  is  the  channel  height  (in  the  y- 
direction)  and  a  the  channel  width  (in  the  z-direction).  The  friction 
coefficient  C£  is  computed  as  a  function  of  the  local  plasma  conditions. 


(3)  Overall  energy  equation 


pu 


d 


h  + 


1  2 
7U 


f 


i  y 


e.n. 

i  i 


(3) 


This  equation  describes  the  transport  of  the  total  energy 


ht  =  h  + 


1  2 
7U 


+ 


7-1 


i=ions 


e.n. 
i  i 


(3’) 


where  e.  is  the  ionization  potential  for  each  ion  i  and  h  the  enthalpy  of 
the  gas  corresponding  to  the  temperature  T.  The  loss  of  energy  due  to 
heat  transfer  to  the  walls,  s^,,  consists  of  two  terms,  i.  e. , 


sE  =  2puSt 


h 


,  -  h 
ad  w 

R 


+ 


V 

s 


IT 


(3H) 
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where  S  is  the  Stanton  number  computed  for  the  local  plasma  conditions, 
*  2 

h^s  h  +  u  where  r  is  the  recovery  factor,  h  the  enthalpy  cor¬ 
responding  to  the  wall  temperature  T  ,  and  V  is  the  sheath  drop  in 
volts.  Note  that  the  sheath  drop  V  constitutes  only  a  part  of  the  total 

S 

electrode  drop  V  j,  the  latter  including  also  the  voltage  drop  due  to  the 
boundary  layers. 


(4)  Electron  energy  equation 

r  d  3  Pe  .  d  1 

pu[_H£7  T  PeH£  p  J  = 

=  B'.  Je  -lknev  6  (Te-Tn)-  Sc.s.-R  (4) 

1 

where  pg  is  the  electron  partial  pressure  pg  =  nekTg  (k  Boltzmann’s 
constant),  E'  =  E  +  UXB,  J  the  electron  current,  v  the  total  col- 

c  t 

lision  frequency  for  electrons  (as  defined  in  Ref.  6),  6e^  the  effective 

energy  loss  factor  for  electrons  (see  Refs.  1  and  2,  and  also  Ref.  9), 

Tn  the  excitation  temperature  of  the  gas  (see  Ref.  10),  and  R  the  radi¬ 
ation  loss  from  the  gas.  Detailed  descriptions  of  each  term  on  the  right- 
hand  side  of  Eq.  (4)  have  been  given  in  previous  work  (Refs.  1-6  and  9- 10). 
Therefore  it  is  necessary  only  to  note  here  that  the  first  term  expresses  the 
Ohmic  heating  of  the  electrons,  the  second  term  the  collisional  energy 
transfer  from  electrons  to  heavy  particles  by  all  elastic  and  inelastic 
collisional  processes  other  than  ionization  or  electronic  excitation,  and 
the  third  term  expresses  the  energy  needed  for  ionization  or  released  by 
recombination. 


Equations  (1),  (2),  (3),  and  (4)  are  the  N  +  3  differential  equations 
used  for  the  solution  of  the  quasi-one -dimensional  problem.  Note  that 
from  Eqs.  (i) 


d 

pU  d£ 


e.n. 


1  1 


s. 

l 


d") 


and  consequently  the  overall  energy  equation  (3)  can  be  written  in  the  form 


d 

pu  7S 


,  ,  i  2 

h  +7U 


=  E 


-  7>.s.  - 
y  i  l 


(3a) 
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The  relation  (1”)  was  used  in  deriving  Eq.  (4). 

In  addition  to  the  above  differential  equations,  two  constitutive 
relations,  Ohm’s  law  and  an  equation  of  state,  are  used  to  express  the 
electric  field  E  and  the  gas  pressure  p,  respectively.  Ohm’s  law  has 
the  following  form  (Refs,  1,  2,  and  6),  in  which  the  contribution  of 
thermal  diffusion  has  been  omitted  for  the  purpose  of  treating  the  core  of 
the  flow: 

E’  =  IT  +  ^T  X  "k  (5) 

where  k  is  the  unit  vector  in  the  z -direction,  and  the  coefficients  tr,|3,e 
are  computed  to  the  second  approximation  following  the  methods  of  Ref.  6. 
On  the  other  hand,  the  equation  of  state  has  the  general  form 

P  =  P<V  Ta)  (6) 

where  n^  and  Tq  are  the  number  density  and  characteristic  temperature, 
respectively,  of  each  component  a.  In  the  present  solution,  the  equation 
of  state  is  used  in  its  simplest  form 

P  =k  neTe  +  (  Yj  nK)T  (6a) 

L  K^e  J 

Finally,  note  that  use  of  the  complete  set  of  Eqs.  (1)  implies  that  the 
overall  conservation  of  mass 

puA  =  m  =  constant  '  (im) 

is  also  satisfied,  where  m  is  the  mass  flow  rate  in  the  channel. 

The  quasi-one -dimensional  formulation  of  the  gasdynamic  problem 
in  the  core  of  the  flow  is  completed  by  specifying  initial  conditions  for  the 
unknowns  at  a  given  station  xq  of  the  channel. 
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2.  2.  Compressible  Turbulent  Boundary  Layer  Development 

Boundary  layer  development  on  the  electrode  walls  of  magneto  - 
hydrodynamic  channels  is  a  major  factor  in  determining  their  performance 
The  serious  problem  of  channel  choking  by  boundary  layer  growth  has 
often  occurred  and  been  reported  in  the  past  (Refs.  11  and  12).  Also, 
Phases  I  and  II  of  this  study  have  demonstrated  that  gasdynamic  boundary 
layers  have  a  decisive  influence  on  the  electrical  performance. 

A  realistic  analytical  treatment  of  these  boundary  layers  must 
account  for  the  fact  that  they  are  compressible  and  turbulent,  and  that 
in  addition  to  the  presence  of  both  mean  and  fluctuating  electromagnetic 
forces  and  heat  sources,  the  very  structure  of  turbulence  is  affected  by 
the  presence  of  strong  electromagnetic  fields  in  these  devices. 

Such  a  treatment  has  been  carried  out  by  this  study,  employing  a 
novel  formulation 'developed  previously  (Ref.  7).  This  formulation 
is  based  on  the  close  relation  that  is  found  t©  exist  between  the  tur¬ 
bulent  shear  stress  t  and  the  turbulence  kinetic  energy.  Using  experi¬ 
mental  data  to  derive  universal  turbulence  structure  parameters  a^, 
and  a^  that  correlate  the  turbulent  shear  stress  with  the  intensity,  dissi¬ 
pation,  and  diffusion  of  turbulence  energy  (see  Ref.  7),  the  turbulence 
energy  equation  (which  is  an  equation  for  the  rate  of  change  of  turbulence 
kinetic  energy  along  a  ..streamline  of  the  mean  flow)  is  converted  into  a 
transport  equation  for  the  rate  of  change  of  turbulent  shear  stress.  Then 
this  equation,  together  with  the  mean  momentum  equation,  the  mean 
continuity  equation  and,  for  the  compressible  case,  the  mean  energy 
equation  form  a  closed  system  (provided  the  turbulent  heat  conduction  q 
can  also  be  related  to  the  other  variables)  that  can  be  solved  numerically. 
Thus,  this  formulation  avoids  the  use  of  phenomenological  "mixing 
length*  expressions,  and  can  treat  the  influence  of  electromagnetic  fields 
on  the  turbulence  structure  quite  naturally:  turbulence  suppression  by 
the  magnetic  field  is  accounted  for  by  an  additional  term  in  the  turbulence 
kinetic  energy  equation. 
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This  novel  formulation  has  been  used  in  this  study  in  parallel  with 
the  traditional  "eddy  viscosity M  approach  so  as  to  provide  a  comparison  of 
the  two,  and  continuous  check  of  the  novel  method.  The  method  of  solution 
for  the  two  formulations  is  basically  the  same  (see  Section  III). 

The  unknowns  of  the  two-dimensional  compressible  turbulent  mag¬ 
netohydrodynamic  bound  ary -layer  problem,  as  formulated  in  this  study, 
are 

(1)  the  mean  velocity  U  =(u,  v), 

(2)  the  turbulent  shear  stress  t  =  -  p(u'v’)  , 

(3)  the  gas  enthalpy  h  =  e  +  p/p, 

(4)  the  number  densities  na  (or  partial  mass  densities  pQ,  or 
concentrations  cq  =  pQ/p)  of  each  component  a,  and 

(5)  the  characteristic  temperatures  Tq  (or  enthalpies  hQ)  of 
each  component  a. 

Primed  symbols  denote  the  fluctuating  part  of  each  quantity;  the  reader  is 
referred  to  Ref.  7  for  detailed  definitions. 

The  formulation  described  in  Ref.  7  presents  a  closed  system  of 
equations  for  the  first  four  of  the  above  unknowns  plus  the  turbulent  energy 
flux  q  s  p(e'v*)  and  the  mass  diffusion  fluxes  d^s  (  pW’)  ,  with  the  require¬ 
ment  that  the  values  of  nine  universal  structure  parameters  to  a^  be  known 
as  functions  of  the  normalized  cross -coordinate  y/8.  However,  sufficient 
experimental  data  are  available  only  for  a^ ,  a^,  and  a^>  only  very  scare* 
measurements  exist  for  the  turbulence  correlations  necessary  for  defining 
the  other  six  structure  parameters.  Consequently,  in  this  study  the  use  of 
turbulent  structure  functions  has  been  limited  to  a^,  a^,  and  a^,  and  for  the 
moment,  the  conventional  concepts  of  turbulent  diffusivity  and  turbulent  heat 
conductivity  have  been  relied  upon,  through  the  use  of  a  turbulent  Schmidt 


The  characteristic  temperatures  T  were  not  considered  as  additional 
unknowns  in  Ref.  7  because  thermaf’equilibrium  (TQ  =  T)  was  assumed 
therein.  Also,  the  energy  e  rather  than  the  enthalpy  h  was  used  in 
Ref.  7  to  characterize  the  thermal  state  of  the  gas;  the  enthalpy  h  was 
chosen  here  because  of  the  availability  of  a  convenient  thermal  equation 
of  state  for  air  (Ref.  13)  in  the  form  h  =  h(p,  T). 
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number  and  turbulent  Prandtl  number,  to  replace  a^  to  a^.  In  addition, 
only  one  characteristic  temperature  other  than  the  static  gas  temperature 
T  has  been  distinguished,  namely  the  electron  temperature  T  . 

The  system  of  equations  corresponding  to  these  unknowns  in  the 
body-oriented  coordinate  system  (x,  y)  are  the  following  (see  Ref.  7): 

(1)  Overall  mas 8  conservation  equation 

4pU  +  TSy  PV  +  ^F<PV>  =  0  W 

The  last  term,  denoting  mass  diffusion  due  to  fluctuations  in  density  and 
transverse  velocity  component,  is  negligible  compared  to  the  other  two 
terms  for  free -stream  Mach  number  up  to  about  5  (because  the  local  Mach 
number,  based  on  the  fluctuating  velocity,  will  then  be  small).  Thus  the 
overall  conservation  equation  will  be  used  in  this  study  in  the  form 

35  pU  +  Vy  pv  =  0  (7a) 

(2)  Overall  momentum  equation 

p{u!;  +  vw}u+<pV>ts7'  w,Ti-  +  T>  = 

*  -  £  +  JyB  <8> 

where  t  is  the  turbulent  shear  stress  and  the  laminar  stress  due  to 
molecular  viscosity  (r^  =  8u/9y).  Note  that  t  »  in  the  main 

part  of  the  turbulent  boundary  layer;  however,  since  tl  is  of  significance 
in  the  wall  region  it  will  be  retained  in  the  momentum  equation.  On  the 
other  hand,  the  term  (p’v1)  ^  is  negligible  compared  to  the  other  terms 
of  Eq,  (8),  for  free -stream  Mach  number  up  to  5,  and  will  be  omitted. 
Thus  the  overall  momentum  equation  will  be  used  in  this  study  in  the  form 

•,{u4  +  vw}u=w[T  +  ,1x.l7]-S  +  Vs  ,8a> 


-15- 


AEDC-TR-70-86 


(3)  Turbulent  shear  transport  equation 
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As  discussed  in  Ref.  7,  the  first  two  terms  on  the  right-hand  side  of 
Eq.  (9)  describe  the  production  and  dissipation,  respectively,  of  turbu¬ 
lence  kinetic  energy  due  to  fluctuations,  while  the  third  term  accounts 
for  turbulence  suppression  by  magnetic  fields,  which  is  due  to  the  inter¬ 
action  between  the  flow  field  and  the  electromagnetic  interactions.  The 
conductivity  cr,  the  Hall  coefficient  (3  and  the  ion-slip  parameter  e  . 
are  defined  in  Refs.  1,  2,  and  6,  while  the  turbulence  structure 
parameters  a^f  a^,  and  are  defined  in  Ref.  7.  Their  values  across 
the  boundary  layer  are  given  in  Section  III. 


(4)  Overall  energy  equation 

+  <p'0  f| +  {p <hV>  +,JL}  * 


dp 
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+  *  =  +  £s(?) 
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(10) 


where  qT  is  the  thermal  conduction  due  to  molecular  heat  conductivity 

J-t 

and  species  diffusion.  The  terms  9qL/9y  and  rL(9u/9y)  are  generally 
small  compared  to  the  other  terms  of  Eq.  (10),  but  will  be  retained  because 
of  their  significance  close  to  the  wall.  On  the  other  hand,  the  term 
(p'v*)  (9h/9y)  is  negligible  for  free -stream  Mach  number  up  to  5  and  will 
be  omitted.  Note  that 


9u  „  /  9u  \ 

TL  77  "  \7y  / 


and 


kT  „  9h  _  9c 

qL  =  '  T~  2  ca  -ST-  '  pDL  S  \  Sy- 

p  a  7  n  7 


(11) 


where  k^,  and  are  the  molecular  viscosity,  heat  conductivity. 
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and  binary  diffusion  coefficient,  respectively.  Using  the  definition 


h  =  Jch 

n  I 


a  a 


(12) 


it  is  possible  to  express 

Z 


e  ®ha  .  8h  y.  8ca 

a  Ca-3T  '  3y  '  “  h»  "W 


(13) 


The  concepts  of  the  laminar  Prandtl  and  Schmidt  numbers  Pr^  and  Sc^, 
respectively,  can  be  introduced  so  that 


c  UT 

pDL  =  3S7 


(14) 


Eqs.  (13)  and  (14)  can  then  be  used  to  express 

qL  =  •  W+  (set  "  r^-)  ] 


(15) 


The  turbulent  heat  flux  p(h,vl)  is  also  expressed  by  using  the 
definition  (12),  which  leads  to  the  result 

P<hV)  =p{E  cQ<hiv’>  +  E  h0<cav'>} 

v.  a  a  J 

=  I7  +  (sk  •  K>)?Vsr]  tl6> 


where  the  concepts  of  the  turbulent  Prandtl  and  Schmidt  numbers  Pr^. 
and  Scrj,,  respectively,  have  been  introduced  to  relate  the  turbulent 
fluxes  p(Vv*)  and  p(c^v')  to  the  gradients  8ha/8y  and  8ca/8y  as 
follows : 

8hQ 

P  (hivI>  =  “  PrTT8u/8y 

r  8cu 

p<civ,>  =  -  gCrrT8u/ffy 


(17) 
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The  turbulent  Prandtl  number  has  generally  been  found  to  have  values 
close  to  0.  9  (Refs.  14  and  15), 

Using  the  expressions  for  q^  and  p^hV*)  given  by  Eqs.  (15) 
and  (16),  respectively,  Eq.  (10)  can  be  written  in  the  form 


Note:  An  alternative  form  of  the  overall  energy  equation,  that  can  be  con¬ 
venient  to  use,  is  the  transport  equation  for  the  total  enthalpy  H  =  h  +  . 

It  has  the  form 

P  {u  4e  +  v^f}H+<PV>  I7  +  w{P<HV>  +<*l}  = 

=  -57  tLu+S''  ^  <18) 

The  turbulent  flux  p(H’v’)  can  be  written 

p^H'v1)  «  p^h'v’)  +  pu(u'v’) 

=  p(h'v1)  -  ut  (19) 

and  again  the  term  (p'v1)  8H/3y  can  be  omitted.  Consequently,  Eq. 

(18)  can  be  written  in  the  form 
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+{|iL(sS£'T4^)+?!i73y  (^  '  P^r)}2ha-37  ]+3’-E’  <18a> 

(4)  Species  conservation  equation 


or,  equivalently, 

f  {“  4  +  v  Sy } ca  +  ( w  -5 7  +  -8?{t‘«v'>  +do,L}  = 

=  ma<  SQ>  (20) 

where  ( sq)  is  the  mean  source  of  particles  of  component  a  due  to 

chemical  reactions.  The  molecular  diffusion  term  3d  T  /3y  is 

a,  J 

generally  small  compared  to  the  other  terms  of  Eq.  (20)  in  the  main  part 
of  the  turbulent  boundary  layer,  but  will  be  retained  because  of  its  sig¬ 
nificance  close  to  the  wall.  It  will  be  expressed  in  terms  of  the  laminar 
Schmidt  number  ScL  as  usual: 

9c  H-  3c 

da,L  =  “  PDL  "Sf  =  ^  (21) 

On  the  other  hand,  the  term  (prv*)  3c  /3y  that  appears  in  the  second 
version  of  Eq.  (20)  is  negligible  for  the  same  reasons  mentioned  in  the 
derivation  of  the  continuity,  momentum,  and  energy  equations,  and  will 
be  omitted. 

The  turbulent  diffusion  flux  p  (  c^v1)  can  be  written  in  terms  of 
the  turbulent  Schmidt  number  Sc^,  as  indicated  by  the  second  of  Eqs.  (17). 

Consequently,  the  second  of  Eqs.  (20)  can  be  given  in  the  follow¬ 
ing  form,  which  will  be  used  in  this  study: 

p{u4z  +  v  =  w(s^-  +  3^T  ^7Sy)"Sy  +mJBJ  (20a) 
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(5)  Electron  energy  equation 

? {“  -L  +  T  ~Sy\  cehe  +  TSy  cehe  + 

+  w{P<(Cehe,,T'>  +  l}  “  {u  4+v^}Pe  = 

=  T  ■  (IT  +  UXB)  +  {  J*  •  (E*  +  U'XB)}  - 

6  6 

-Ikn  v.  6  .,(T  -T  ).  J  £.(•.)  -  R 
c  e  t  eff'  e  n'  .  ^  i'  v 

ions 


(22) 


where  cg  =  pe/p  is  the  electron  concentration,  hg  =  (5k/2me)Te  the 

electron  enthalpy  (c  =  5k/ 2m  },  p  =  n  kT  the  electron  pressure, 

p,  e  e  e  e  e 

Te  the  electron  temperature,  the  excitation  temperature  of  the  gas 

(Ref.  10),  the  ionization  potential  for  each  ionizing  species,  and 
q  T  the  heat  flux  due  to  conduction  and  diffusion  of  electrons,  expressed 

e,  -Li 

as 


3T 


0c 


le,h  =  "  ke  “5y  "  S? 


fce  9he  . 

=  - 


1  h  8ce  I 

szr\-&Tj 


(23) 


where  k^  is  the  thermal  conductivity  for  electrons.  Note  that  Eq.  (23) 
is  consistent  with  the  multicomponent  expression  for  qL  given  by 
Eq.  (11).  Again,  for  free-stream  Mach  number  up  to  5,  the  term 

cehe  is  negligible  compared  to  the  other  terms  of  Eq.  (22)  and 
will  be  omitted. 

The  turbulence  flux  p((c  h  J'v1)  can  be  expressed  as  follows 

©  © 


p<(c  h  )V>  =  pc  (h* v*)  +ph  (c'v*) 


e'  e 
0u/dy 


e'  e 
8h 


©  C  i,  C  G 

VxZ  &y  ScT  cJy 


(24) 


-20- 


AEDC-TR -70-86 


where  use  has  been  made  of  Eqs,  ( 17). 

Finally,  the  total  current  density  J  in  a  magnetohydrodynamic 
channel  is  very  closely  approximated  by  the  electronic  current  J  alone, 
so  that  the  results  of  Ref.  7  can  be  used  to  express 


(  •  (E1  +U'XB))«  (J'.  (E1  +  U'  X  B  )>  = 


1  eo-BZ  t 
3  e2+|3*  aiP 


(25) 


Using  the  expressions  given  by  Eqs.  (23)  -  (25),  Eq.  (22)  can  be 
written  in  the  form 


Piu4  +  v|F}cehe  = 


9  [  /  x  T  1  \  9he  ,/^L  +  T  1  V  9cel  + 

szi&f  vrT)  ce  +\ 3^;  +  -511757  + 

+  {“4  +  v|fK  +7e’  <E+ffxB)+  3 


e<r  B  t_ 

T  eZ+pZ  alP 


(22a) 


10ns 


Equations  (7a),  (8a),  (9),  (10a)  or  (18a),  (20a)  and  (22a)  constitute 
a  system  of  N  +  5  partial  differential  equations  for  ttye  N  +  5  unknown 
quantities  u,  v,  t,  T,  c  and  T  of  the  gasdynamic  problem  in  the 
compressible  turbulent  boundary  layers. 

The  use  of  a  transport  equation  for  t  to  replace  the  conventional 
algebraic  "eddy-viscosity"  relation  is  discussed  in  detail  in  Ref.  7.  It 
results  in  a  momentum  equation  [Eq.  (8a)]  in  which  the  diffusion  term 
3r/3y  is  not  of  gradient  diffusion  type.  It  is  worthy  of  note  that  the 
diffusion  term  in  the  turbulent  shear  stress  equation  (9)  is  not  of  gradient- 
diffusion  type  either  (Ref,  7).  As  a  consequence  of  these  two  facts,  the 
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system  of  governing  equations  is  hyperbolic. 

Note  that,  when  convection,  diffusion,  and  magnetic  suppression 
terms  are  negligible,  Eq.  (9)  reduces  to  an  algebraic  expression  of  the 
form  of  the  mixing  length  relation,  namely 


(t/p) 


1/2 


a26 


8u 

*7 


(9a) 


with  the  dissipation  length  appearing  in  place  of  Prandtl's  mixing 

length  f ,  This  situation  prevails  close  to  the  wall.  On  the  other  hand, 
near  the  outer  edge  of  the  boundary  layer  the  diffusion  and  convection 
terms  are  dominant  in  Eq.  (9),  which  therefore  assumes  the  form 


{u  4  * v 


zi^p 


(V  p) 


max 


IT 


00 


a 

■57 


a3T=  0 


(9b) 


This  form  of  Eq.  (9)  is  utilized  to  relate  the  entrainment  rate  m^, 
through  the  outer  boundary  to  the  turbulence  structure  parameters 
and  Let  the  outer  edge  of  the  boundary  layer  be  defined  as  the  line 
yE(x)  on  which  the  streamwise  velocity  u  equals  99.5%  of  the  free - 
stream  value;  in  other  words,  yE(x)  =  6g^g(x),  where  u(x,  6995)  = 

=  0.  995  at  every  x.  Accordingly,  Eq.  (9b)  yields 

m—  =  -  2a,  a- p  U  "  *  (  — \  (i 

E  1  3roo  oo  \  p/ 

v  max 


Eq.  (26)  has  been  used  to  obtain  the  structure  parameter  a^(y)  experi¬ 
mentally  through  measurements  of  the  entrainment  rate  (Ref.  14). 
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HI.  METHODS  OF  SOLUTION 

1.  ELECTRICAL  PROBLEM 

The  solution  of  the  electrical  problem  in  the  entrance  and  exit 
regions  follows  closely  the  methods  used  in  Phases  I  and  n  for  the  single 
electrode  problem.  In  fact,  the  solution  of  the  streamfunction  equation 
employs  exactly  the  same  method,  with  different  boundary  conditions, 
while  the  simultaneous  solution  'of  the  continuity  and  electron  energy 
equations  on  each  streamline  (see  Ref.  2)  is  simplified  by  the  fact  that 
initial,  rather  than  boundary -type,  conditions  are  now  used  for  nQ  and 
T  .  (This  fact  eliminates  the  need  for  carrying  out  the  "initial  condition 
iteration"  described  in  Ref.  2,  ) 

The  rectangular  variable -me  sh  grid  of  points  used  for  the  numeri¬ 
cal  solution  consists,  as  in  the  single  electrode  problem,  of  m  rows 
and  n  columns,  where  row  i  coincides  with  the  anode  surface  y  =  0 
and  row  m  coincides  with  the  cathode  surface  y  =  D.  However, 
columns  1  and  n  now  define  the  upstream  and  downstream  boundaries 
of  the  whole  entrance  (or  exit)  region.  Care  is  taken  that  each  of  the 
points  on  the  y-walls  that  separate  conductor  from  insulator  segments 
(i.  e. ,  each  of  the  singular  points  of  the  geometry)  lies  on  a  grid  column; 
this  is  always  possible  owing  to  the  variable  spacing  of  the  grid. 

As  mentioned  above,  the  finite -difference  solution  of  the  stream- 
function  equation  in  the  entrance  and  exit  regions  is  accomplished  by  the 
same  direct  method  that  was  employed  in  Phases  I  and  II  (refs,  i  and  2). 
The  number  of  "fundamental  unknowns"  (see  Ref.  1)  is  now  n  -  2. 

Finally,  concerning  the  simultaneous  solution  of  the  continuity 
and  electron  energy  equations,  it  is  well  to  elaborate  on  .the  important 
point  mentioned  above:  When  solving  the  electrical  problem  in  the 
periodic  part  of  the  channel  without  any  prior  knowledge  of  the  plasma 
properties  at  the  initial  plane  of  one  electrode -pair  region  (as  ^»s  done 
in  Phases  I  and  II),  boundary-ttype  (i.  e. ,  periodicity)  conditions  with 
with  respect  to  x  must  be  used,  both  for  the.  current  streamfunction  and 
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for  the  plasma  properties.  But  in  solving  the  electrical  problem  for  the 
entrance  and  exit  regions,  initial  type  conditions  for  nQ  and  T  ,  and 
for  the  plasma  properties  can  be  used.  Furthermore,  since  after 
solution  of  the  entrance  problem  the  values  of  nQ  and  Tg  will  have 
been  obtained  at  the  beginning  of  the  periodic  part  of  the  channel,  initial 
type  conditions  for  nQ  and  Tg  can  be  used  throughout  the  downstream 
solution  of  the  electrical  problem  for  subsequent  electrode  pairs.  This 
leads  to  a  significant  reduction  in  computation  time. 

Appendix  A  describes  the  computer  program  INLET,  which  has 
been  coded  in  FORTRAN  for  the  solution  of  the  electrical  problem  in 
the  entrance  or  exit  regions.  It  includes  a  list  of  routines  and  important 
variables,  specification  of  required  input  data  and  formats,  and  a 
description  of  available  output. 

2.  GASDYNAMIC  PROBLEM 

In  solving  the  two-dimensional  gasdynamic  problem  the  flow 
region  has  been  separated,  as  mentioned  in  Section  II,  into  the  core  and 
the  boundary  layers,  at  least  up  to  the  point  where  the  latter  may  meet 
to  result  in  fully  developed  viscous  flow.  The  gasdynamic  problems 
of  the  core  flow  and  the  boundary  layers  are  coupled  (in  addition  to 
their  being  coupled  to  the  results  of  the  electrical  problem):  the  boundary 
layer  growth  determines  the  shape  of  the  effective  channel  area  A(x)  to 
be  used  in  the  core  flow  solution,  while  the  impressed  free  stream  con¬ 
ditions  determine  the  boundary  layer  structure  and  growth  rate. 

2.  1.  Core  Flow  Solution 

The  solution  of  the  quasi-one -dimensional  problem  for  the  core 
flow,  which  was  formulated  in  Section  II,  is  accomplished  by  a  straight¬ 
forward  4th  order  Runge-Kutta-Gill  integration  method.  In  this  formu¬ 
lation  the  magnetic  induction  B  in  the  z -direction,  the  component  J 
of  the  current  density  in  the  y- direction,  and  the  effective  channel 
area  A  are  specified  functions  of  the  streamwise  coordinate  x. 
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For  the  purposes  of  this  solution,  the  ordinary  differential  equa¬ 
tions  of  the  problem,  given  in  Section  II,  are  put  in  canonical  form, 
namely 

"y’  =  T(x,  ‘y  )  (27) 

where  the  vector  y  (x)  has  as  components  the  N  +  3  unknowns  of  the 
problem  at  each  station  x,  while  y 1  is  the  vector  ,  and  7  the 

appropriate  vector  function.  The  vector  function  7  is  readily  found 
by  transformation  of  Eqs.  (1)  -  (4)  of  Section  II.  Employing  also  the 
overall  continuity  equation  { 1 m ) 


puA  =  m  =  constant. 


the  components  of  f  are  then  given  by  the  following  expressions 
(1)  Continuity  equations  for  each  component  a 
s 


d 

Hx 


a  /  1  dA  ,  i  du\ 
n  =  —  -  n  i  ~t  j  *  —  ■t—  | 
a  u  a\  Adx  u  dx/ 


(28) 


(2) 


(3) 


Overall  momentum  equation 


du 


\RM 

u(M2-i) 


(S,  +sz) 


Overall  energy  equation 


dT 

dx 


(yM2-l)S1  +  (y-i)M2S2 
M2-  1 


(29) 


(30) 


(4)  Electron  energy  equation 


d  T  . 

dx  e 


(2/3k)(E’-  7e-  Je.s.)  -  Ve 


n  u 

e 

vt6eff<VTn>  x2 

u  3 


<3‘> 


where 
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ir-r-sE+(VW/D)Jy-^Ciei 

1  = 

SM  '  JyB  T  dA 

2  -  - pi - X 


(32) 

(33) 


and  R  is  the  gas  constant,  v  =  c  /c  the  ratio  of  specific  heats,  M 

,1/2  P  v 

=  u/  (yRT)  the  Mach  number,  and  the  part  of  the  electrode 

voltage  drop  Ve^  that  is  due  to  the  boundary  layers.  All  other  quan¬ 
tities  have  been  defined  in  Section  II. 


The  Runge-Kutta-Gill  integration  formula,  giving  the  value  of 
y  at  the  point  when  its  value  at  the  point  has  been  found, 

reads 


^n+i 


=  y. 


n 


1 

' E 


+  2(1  i/2  yk2  +  2(i+yji/2y£3  +  kz 


(34) 


where  yn+1  =  y  (xn+1),  Yn  =  Y  (xj,  *-n+l  =  xn  +  h  fh  the  integration 
step)  and 

k,  =  h  f  (x  ,  y  ) 

1  '  n’  7  n' 


h  f  (x 


n 


+  ih. 


rn 


1 

7 


kt) 


(34a) 


hT(xn+ih,7n+  (-i  +  >/i72)k1  +(i  -yji/2)kz) 


~k4  =  h  7 (xn+h,  7n  -  yfUI  "k2  +  (1  +  yJT/2)k3  ) 


5 

The  error  in  this  integration  formula  is  of  order  0(h  ).  In  practice, 
accuracy  is  controlled  by  varying  h  so  as  to  limit  the  relative  changes 
of  all  unknown  quantities  over  each  integration  step  to  given  maximum 
values. 


The  above  system  of  differential  equations  is  a  "stiff"  set  of 
equations,  for  the  low  mass  of  the  electrons  results  in  a  comparatively 
small  time  constant  for  Eq.  (31).  Because  of  this  disparity  in  time 
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constants,  it  is  often  desirable  to  treat  the  electron  energy  equation 
separately,  by  replacing  the  x-derivative  with  a  two-point  difference 
expression  and  reformulating  Eq.(31)  into  an  implicit  algebraic  equation, 
which  may  be  solved  by  iteration  (see  Ref.  2).  The  resulting  equation 
is 


Te<X>  = 


J*  -  A  €.8. 
e  y  1  1  „  y 

1  +  3  ( 

-A  T1  4-  U  T-  \ 

n  k  7  \ 

t  eff  x  -  x  e,  of 

e 

o  * 

3  (v  6  +  u  +Se> 

+  du  + 

u 

dA 

IT  \vt° eff  x-x  n  ; 

o  e 

'  dx 

3. 

(31a) 


where  xq  is  the  value  of  x  at  the  beginning  of  the  current  integration 
step,  and  Teo  =  Te(xo). 

Appendix  B  describes  the  computer  program  CHANNEL,  which 
has  been  coded  in  FORTRAN  for  the  solution  of  the  gasdynamic  problem 
in  the  core  region.  It  includes  a  list  of  routines  and  important  variables, 
specification  of  required  input  data  and  formats,  and  a  description  of 
available  output. 


2,  2.  Boundary  Layer  Solution 

For  the  solution  of  the  compressible  turbulent  magnetohydrody¬ 
namic  boundary  layer  problem,  which  was  formulated  in  Section  II,  a 
finite-difference  method  has  been  developed  incorporating  some  of  the 
desirable  features  of  the  method  originally  suggested  by  Patankar 
(Ref.  15)  for  the  treatment  of  the  incompressible  eddy- viscosity  model. 

As  developed  in  this  work,  the  method  of  solution  is  applicable 
to  compressible  or  incompressible,  two-dimensional  or  axi-symmetric 
flows  over  straight  or  slightly  curved  surfaces,  with  or  without  mag¬ 
netohydrodynamic  effects,  and  with  or  without  mass  injection  or  suction 
at  the  walls.  The  most  attractive  characteristics  of  this  method  are: 

(1)  The  use  of  the  normalized  streamf unction  as  the  cross - 
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coordinate.  This  results  in  a  coordinate  system  that  adjusts  automati¬ 
cally  to  the  thickness  of  the  boundary  layer,  while  the  number  of  calcu¬ 
lation  steps  in  the  transverse  direction  remains  constant  as  the  boundary 
layer  grows.  In  other  words,  the  step  size  in  the  mathematical  coordi¬ 
nate  system  remains  constant,  although  the  physical  distance  corre¬ 
sponding  to  this  step  size  varies  as  the  boundary  layer  grows  in  the 
streamwise  direction.  This  is  in  marked  contrast  with  the  complicated 
artificial  devices  that  are  usually  employed  to  satisfy  the  conflicting 
requirements  of  covering  the  entire  boundary  layer  without  allowing 
the  number  of  transverse  calculation  steps  to  grow  excessively  with 
growing  boundary  layer  thickness  (Ref.  18). 

(2)  The  use  of  an  implicit  difference  scheme  for  the  numerical 
solution  of  the  boundary  layer  equations.  Use  of  this  scheme  assures 
stability  of  the  numerical  method  irrespective  of  step  size.  In  other 
words,  step  size  can  be  chosen  so  as  to  satisfy  a  given  requirement 

of  numerical  accuracy  without  concern  about  the  stability  of  the  method. 

(3)  The  use  of  universal,  Couette-type,  solutions  for  the  thin 
laminar  sublayer  (adjacent  to  the  wall).  In  this  region  of  extremely 
steep  transverse  gradients,  use  of  a  finite -difference  computation  would 
necessitate  extremely  small  step  sizes  and  thus  excessive  computer 
time.  Fortunately,  in  this  region  one  can  take  advantage  of  the  fact 
that  inertia  forces  are  very  small  compared  to  viscous  stresses  and 
pressure  forces,  so  that  the  boundary  layer  equations  can  be  given  the 
form  of  total  differential  equations  with  derivatives  only  in  the  cross - 
stream  directions.  This  feature  is  used  in  the  present  method  and  leads 
to  significant  reduction  in  computer  time  compared  to  previous  methods. 

2.  2.  t.  Coordinate  system 

The  coordinate  system  used  is  shown  on  Fig.  3.  It  can  be  seen  that 
both  axi-symmetric  and  plane  flows  are  allowed  for,  with  x  in  the  approxi¬ 
mately  streamwise  direction  and  y  in  the  transverse  direction.  The 
transverse  geometrical  coordinate  y  is  measured  from  the  I-boundary  of 
the  layer.  The  direction  of  the  I-streamlines  (i.  e. ,  the  direction  of  the 
inner  wall)  is  allowed  to  vary  with  x;  however,  this  variation  is  assumed 
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Fig.  3.  The  coordinate  system. 
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to  be  sufficiently  slow  so  as  not  to  give  rise  to  any  additional  (centrifugal, 
or  Coriolis)  forces. 


The  flow  is  assumed  to  be  bounded  by  the  I  and  the  E  boundaries. 
The  discussion  here  will  concentrate  on  the  case  where  the  I-boundary 
is  a  flat  material  wall,  and  the  E -boundary  the  edge  of  the  boundary 
layer,  i.  e.  a  free  boundary.  The  method  can  easily  be  adapted  to  other 
flow  situations  (e.  g.  mixing  layers). 


The  coordinate  system  is  transformed  by  using  the  von  Mises 
transformation  (Ref.  19).  This  is  done  by  introducing  the  streamf unction 
S,  defined  by 

pU  =  VX  (S"k )  (35) 


where  again  k  is  the  unit  vector  in  the  z -direction.  Eq.  (35)  of  course 
implies  that 


as 

W 


=  +  pu 


(35a) 


so  that  the  continuity  equation  (7a)  is  satisfied  automatically.  All  equations 
cin  notv be- transformed  -from  the  system  (x,  y)  to  the -system  (x,  S)  by 
it's  mg  the  relations 


(s)y-  (^)s  ■ pv  (4)a 

(w)x = pu(^)x 

Under  this  transformation,  the  convective  operator 
the  (x,  y)  coordinate  system  has  the  form  p 
becomes  simply  p  u(^:)  »  which  is  of  course' 
definition  of  the  streamfunction. 


(36) 


pU  •  V  which  in 

P  ax'! 

now 
the 


erator  pU  •  V  which 

{u(X)y +  v(^)} 

e  a  direct  result  of  t) 


Following  the  von  Mises  transformation,  the  coordinate.  S  is 
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normalized  by  again  transforming  the  coordinate  system  from  (x,  S)  to 
(x,co)  where 


u 


S  -  S, 


SE  "  SI 


(37) 


and  S^(x)  and  Sg(x)  are  the  values  of  the  streamf unction  at  the  two 
boundaries  of  the  flow  for  a  given  x.  As  long  as  the  possibility  of 
mean  flow  reversal  is  excluded  (which  is  actually  the  condition  for  the 
streamfunction  to  be  single -valued),  the  boundary  layer  flow  region  is 
always  imbedded  in  the  constant  interval  0  <  cj  <  1. 


It  is  clear  from  the  definition  of  Sj  and  SE  that 

^  SI  =  '  mI 

d  q  _  • 
cGc  SE  ~  ‘  mE 


(38) 


where  m^  and  mE  are  the  mass  transfer  rates  across  the  two  bound¬ 
aries.  Consequently,  the  transformation  from  the  system  (x,  y)  to 
the  system  (x,  to)  will  proceed  according  to  the  relations 


-pv 


+ 


m.j  +  w(mE 

BE-Si 


(39) 


Accordingly,  the  convective  operator  now  becomes 


plf-v  =  pu[(4)^  + 

where 

mT  +  u  (m^  -  mT) 

R(x,  w)  =  - £ - L 

SE  "  5I 


(40) 


(40a) 
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2.  2.  2.  The  transformed  equations 


In  the  "modified"  von  Mises  coordinate  system  (x,  w)  the  govern - 
ing  equations  of  the  problem  (see  Section  II)  read  as  follows: 


( 1)  Overall  momentum  equation 


+r  ■L'y  - 


3  _l  „  3  L  3t/3w  i 


E  “I  (SE-  Sj) 


3  .  3u 

2  333  ^P"  333  = 


=  i-{. 

pu  t 


P  +  J  B  f 
dx  y  j 


(2)  Turbulent  shear  transport  equation 

(T/P>max  9 


{353 +  R  ^3}  Wrj 


^SE"  Sj)  333  a3T  = 


T  3u 


-i-fj-fjc's  + 

pu|_a?6  \  p  / 


,3/2  ,  _  „2 


1  eo-B 


2-1 

llP  J 


(3)  Overall  energy  equation 


(41) 


(42) 


/  9  4.  9  \  i,  1  9  f/^L  ,  T  \  9h  , 

X&  +  R  333jh  =  ~  ~-j'2  333  pu|_\  VT^  Ft "  3379y  >  TO 

E  I 

+{|1l(s^;  "t4^)  +  ■5V57  (5^;  •  PF^)  }2ha-5^-]  + 

+  1  &  +  -!££.- (»£)*  +  J_r^_(lf2  + 

o  dx  /e  c,  .2  \  aw/  pu  a_6  \  p  / 


<Vsi> 


pu  j_a26  \  p 

+  7  •  (E+UXB)  +  i  - 

5  €^+(3^  a 


T 

i  P  J 


(43) 


or,  for  the  total  enthalpy  H, 
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TSpsjT 


_u{(1“P?l)T  + 


W*  9u' 
SE‘^I  ’85- 


If.? 

pu 


(4)  Species  conservation  equation 


(5)  Electron  energy  equation 


Equations  (41),  (42),  (43)  or  (44),  (45),  and  (46)  constitute  the 
system  of  partial  differential  equations  that  are  solved  numerically 
for  the  unknowns  u,  t,  T,  cq,  and  Tg.  Note  that  the  density  p  is  not 
a  true  additional  variable,  since  it  is  related  to  the  temperature  and 
pressure,  the  latter  being  a  given  function  of  x  determined  by  the  free- 
stream  conditions. 
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2.  2.  3.  Finite  difference  formulation 

An  important  step  in  the  numerical  method  is  the  discretization 
procedure,  i.  e.  the  transformation  of  the  differential  equations  of  the 
problem  into  a  suitable  system  of  algebraic,  finite -difference  equations. 

To  derive  these  finite -difference  expressions  an  implicit 
scheme  suggested  by  Pakankar  has  been  utilized  (Ref.  17).  The  notation 
appropriate  to  this  derivation  is  shown  on  Fig.  4.  A  rectangular  grid 
of  points  is  used  in  the  x-o>  plane  with  variable  spacing  in  both  the  x- 
and  to- directions,  to  provide  for  desirable  requirements  of  accuracy 
in  regions  of  special  interest,  e.  g.  near  separation,  or  near  an 
electrode  spot.  This  grid  requires  no  adjustment  in  the  <o-direction  in 
succeeding  x-stations.  It  is  assumed  that  the  values  of  all  dependent 
variables  are  known  at  the  grid  points  U",  U,  and  U+,  which  belong  to 
the  upstream  station  x  =  x^.  Their  values  at  the  points  D“,  D  and  D+, 
which  belong  to  the  downstream  station  x  =  x^  are  to  be  found.  The 
finite  difference  equations  are  derived  by  integrating  each  of  the  differ¬ 
ential  equations  of  the  problem  over  a  control  volume  bounded  by  the 
points  UU+,  UU~,  DD+  and  DD“  under  the  assumption  of  a  simple  vari¬ 
ation  of  all  dependent  variables  between  grid  points.  Specifically,  in 
the  w-direction  the  dependent  variables  are  assumed  to  vary  linearly  (or 
quadratically)  between  grid  points,  while  in  the  x  direction  a  step-wise 
(or  linear)  variation  is  assumed.  When  a  step-wise  x- variation  i« 
utilized,  then  the  values  of  the  dependent  variables  in  the  interval 
between  xu  and  x^  (except  at  xu)  are  assumed  uniform  and  equal 
to  those  at  x^.  The  use  of  this  procedure  guarantees  that  the  conser¬ 
vation  equations  will  be  satisfied  uniformly  throughout  the  boundary  layer. 
Note  also  the  use  . of  a  step*wis*<  x-varilatibn  implies  that  tibia  derivatives 
in  the  u-direction  are  evaluated  at  x  =  xd.  This  implicit  difference 
scheme  is  known  to  be  unconditionally  stable  for  linear  systems  of 
differential  equations  (Ref.  20).  Of  course,  the  system  of  equations 
governing  this  problem  is  not  a  linear  one.  Nevertheless,  the  employed 
implicit  difference  scheme  has  shown  good  stability. 

The  important  question  that  arises  when  attempting  to  solve  the 
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Control 

Area 


Fig.  4.  Control  area  for  the  implicit  finite -difference  scheme. 


Fig.  5.  Definition  of  "slip  values'*  for  each  dependent  variable 
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system  of  derived  finite -difference  equations  is  how  to  maintain  a  linear 
system  (which  is  readily  solvable)  without  loss  of  accuracy.  This  is 
accomplished  in  this  method  of  solution  by  combined  application  of  the 
techniques  of  quasi-linearization  and  iteration;  in  other  words,  the  known 
upstream  values  of  the  dependent  variables  are  utilized  to  evaluate 
selected  coefficients  in  the  system  of  equations  so  as  to  linearize  it, 
and  then  the  results  are  enriched  by  iteration. 


The  governing  differential  equations,  except  Eq.  (42)  where  the 
diffusion  term  is  not  of  gradient -diffusion  type,  have  the  general  form 


8$ 


+  (a  +  bo) 


0®  a  0®  .  , 

c)a>  =  C  c>u  d 


(47) 


where  the  left-hand  side  expresses  the  convection  term  and  the  right- 
hand  side  contains  the  diffusion  term  and  a  general  source  term  d. 

It  is  clear  that  the  convection  term  is  linear,  involving  linear  first 
order  differential  operators,  and  therefore  results  in  linear  finite - 
difference  expressions  for  each  unknown  ©.  On  the  other  hand,  the 
diffusion  term  is  not  always  linear  with  respect  to  the  dependent  vari¬ 
ables  because  c  can  be  a  function  of  some  of  them.  Similarly,  the 
source  term  d,  involving  general  production  and  dissipation  terms,  is 
generally  nonlinear  with  respect  to  the  dependent  variables.  The 
treatment  of  the  source  terms  in  the  solution  procedure  is  extremely 
critical  because  of  the  importance  of  production  and  dissipation  terms 
in  the  boundary  layer  equations.  A  method  of  quasi-linearization  has 
been  successfully  adopted,  which  entails  expressing  the  downstream 
value  d^  of  each  source  term  in  terms  of  its  known  upstream  value 

d  as  follows : 
u 

d,  =  d  +  (S’.  - F  )  •  (V,d)  (48) 

d  u  d  u  ©  u 

where  ©  is  the  vector  of  the  unknowns  (u,  t,  T,  ca>  T  ). 

When  all  terms  in  the  differential  equations  of  the  problem  are 
transformed  by  the  quasi-linearization  procedure  described  above,  then 
to  each  control  volume  across  the  boundary  layer  at  a  station  x  there 
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corresponds  a  system  of  linear  finite  differences  equations  of- the  form 

=  A  •  7  .  +  B  .  S’  +C  (49) 

D  D  D' 

where  A  and  B  are  matrices  and  C  a  vector  of  coefficients.  There 
are  n-2  such  control  volumes  across  the  boundary  layer,  where  n  is  the 
number  of  grid  points  in  the  w-direction.  On  the  other  hand,  there  are  n 
unknowns  $  .  The  additional  2  equations  needed  for  the  solution  will  be 
the  boundary  conditions  at  the  two  edges  of  the  boundary  layer.  At  the 
outer  edge  <o  =  i,  the  boundary  condition  may  be  applied  simply  by 
setting  the  values  of  the  dependent  variables  at  that  point  equal  to  their 
free-stream  values.  At  the  inner  edge  w  =  0,  on  the  other  hand,  the 
boundary  condition  will  be  expressed  by  utilizing  Couette-type  solutions, 
which  will  account  for  the  .strong  gradients  of  the-  dependent  variables 
in  the  wall  region. 

2.  2.4.  The  wall-region  Couette-type  solutions 

It  has  already  been  mentioned  that  an  important  feature  of  this 
method,  leading  to  significant  savings  in  computer  time,  is  the  treat¬ 
ment  of  the  wall- region  of  the  boundary  layer,  i.  e.  the  thin  "sublayer" 
adjacent  to  the  wall. 

There  are  two  main  reasons  for  a  special  treatment  of  this 
"sublayer": 

(1)  Since  the  scale  and  the  intensity  of  turbulence  tend  to  zero 
close  to  a  solid  wall,  the  laminar  contributions  to  the  fluxes  become 
significant,  and,  very  close  to  the  wall  (in  the  so-called  "laminar 
sublayer"),  they  become  dominant.  This  applies  both  to  momentum 
transport  and  to  heat  conduction. 

(2)  The  transverse  gradients  of  the  dependent  variables  become 
very  steep  in  this  region,  so  that  a  finite -difference  solution  in  this 
region  would  require  an  excessively  fine  grid  in  the  w-direction. 

This  method  takes  advantage  of  the  fact  that  in  the  wall  region, 
where  u  is  small,  streamwise  convective  effects  (and  thus  stream- 
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wise  derivatives)  can  be  neglected  in  the  governing  equations.  In  other 
words,  there  exists  a  one -dimensional  boundary  layer  near  the  wall; 
this  is  what  is  also  known  as  "Couette  flow.  w  Then  the  flow  problem  in 
this  region  reduces  to  solving  ordinary  rather  than  partial  differential 
equations,  and,  once  these  solutions  are  obtained,  they  will  be  uniformly 
valid  (i.  e.  independently  of  x).  This  means  that  these  "inner  solutions" 
can  then  be  used  as  asymptotic  forms,  or  boundary  conditions,  for  the 
solution  of  the  partial  differential  equations  in  the  main  part  of  the 
boundary  layer. 

In  obtaining  the  Couette-type  solutions  in  the  wall-region,  the 
turbulent  shear  stress  r  is  expressed  by  Eq.  (9a),  which  is  valid  in  this 
region  where  diffusion  and  convection  of  turbulence  energy  are  indeed 
negligible.  Thus,  the  concept  of- the  Peddy- viscosity"  can  be  applied 

in- this  region  without  contradicting  the  transport  formulation  for  the  tur¬ 
bulent-shear,  andEq.  (42)  can  be  replaced  by  the  algebraic  relation 

I^J,  =  p(a.2sj2  j  -gH 

=  W57  ,50) 

Under  these  assumptions,  the  following  Couette -type  -first  integrals 
of  the.. remaining  governing,  equations  can  be  obtained: 

(1)  Overall  momentum 

(ixl  +  M H7  =  Tw+  V  +  y  -  /J  JyB  (51) 

where  m^  is  the  mass  flux  through  the  well  due  to  injection  or  suction, 
the  shear  stress  at  the  wall,  and  the  subscript  w  generally  denotes 
conditions  at  the  solid  wall. 

(2)  Overall  energy  equation 

/^L  ,  ^T  \dh  „  .  •  ,,  ,  .  , 

(P5^  +  Pr^)cly  =V  mw(h"V  + 

+\  “  3^)  + 
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or 


-  K  ■  ^i)  ?i'“'3f}w'  -C  [ u  s + ^(H)  + 

+  S ^r(?)  +  7  •  <E  +U  XB)  +  ^  ]dy 

(P^  +  P^)  If  +  U  {  (‘  ‘  PJ^)^  +  (‘  -  P^,)1^}  37  " 

■  -v  +  mjH- V  + 

+{  •‘i/ps^  ’  +  >*t(p^.  -  5E^)}  l 

-{"iXps:  -  sej)  2  ha  -ar}w  -  T  dy 


dc 

a 

d  dy 


(52) 


(53) 


(3)  Species  continuity  equation 
dc 

c 

'L  “'-ip/  3y 


/  l^T  \  dC  ■  rY 

(■* —  +  •« — )  -3 —  =  d  +m(c  -c  )  -  m  J  ( s  )  dy 
\ScT  ScZ/  dy  a,w  w'  a  a,w'  aJ  '  of  y 


(54) 


(4)  Electron  energy  equation 
dh 

)Co 

‘L  '  ‘T 


/tV  \  dh 

(pj£  +  PE^.  )  ce  -3T  =  qe,w  +  (mwce,w"  de,w>  ’  <V  Vw*  ‘ 

-  /Y[v^  +7e'  &  +  UX5)  + 


TT  V/  Y?  V  _L  1  ^ 

^  €Z+p2 


T 

a{P 


3 

7 


knevt^eff(Te"  Tn^  “  .£  {(ei  "  meheK8i>}  -  R  ]  dy  (55) 

10ns  j  J 


where  the  species  continuity  equation  (54)  for  electrons  (a  s  e)  has  been 
used. 


In  order  to  solve  these  first  integrals,  to  obtain  the-  Couettef- 
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type  solutions  for  the  unknown  vector  <j>  =  (uf  T,  c  ,  T  ),-•  an  analytical 
expression  for  or  a26  is  necessary,  that  will  be  valid  and  accurate 
very  close  to  the  wall.  For  this  purpose  von  Karman's  expression  Ky 
(K  =  0,4  is  von  Karman's  constant),  multiplied  by  an  exponential  damping 
factor,  is  used -for  the  dissipation  length  a26  in  the  wall  region,  as 
suggested  by  Van  Driest  {Refs.  14  and  21): 


a26  =  Ky[  l  exp  j-  3^-  J  V H  V>} 


(56) 


The  values  0.  4  and  26.  0  have  been  employed  for  the  dimensionless 
constants  K  and  At,  respectively  {Ref.  21).  The  presence  of  the 
exponential  factor  in  Eq.  {56)  allows  its  use  arbitrarily  close  to  the 
wall,  i.  e. ,  even  in  a  purely  laminar  region. 

It  is  then  possible  to  integrate  Eqs.  (51)- (55)  once  more  to 
obtain  the  Couette  flow  solutions  for  the  dependent  variables  6.  In 
general,  this  second  integration  will  have  to  be  carried  out  numerically 
to  yield 

7  =?<<■>.  7W.  Jjg  ,  mw,  T,  E,  B)  (57) 

where  the  vector  F  =  (  t  ,  q  ,  d  ,  q  )  contains  the  wall  fluxes  cor- 

w  w’  w  a.  w  e,  w 

responding  to  =  {u,  T,  ca,  Te).  At  a  given  station  x,  the  mass  flux 
through  the  wall  m^  and  the  wall  value  <{>w  of  the  dependent  variables  is 
generally  known.  In  addition  the  pressure  gradient  dp/dx  and  the 
electrical  fields  J,  E  and  B  are  prescribed.  Thus,  the  Couette-type 
"inner"  solutions  will  contain  as  parameters  the  wall  fluxes  F  .  The 
values  of  these  parameters  will  be  determined  by  matching  both  the  mag¬ 
nitude  and  the  derivative  of  the  "inner"  to  the  "outer"  solutions  at  some 
appropriate  points  y^,  where  both  solutions  remain  valid.  Fortunately, 
the  choice  of  y^,  is  not  critical,  since  both  solutions  are  valid  within  a 
fairly  broad  common  domain,  the  so-called  "log"  region  or  "buffer" 
layer  that  lies  between  the  thin  laminar  sublayer  and  the  region  where 
streamwise  convection  effects  start  being  significant.  In  fact,  due  to 
the  extreme  thinness  of  the  inner  region,  it  is  usually  sufficient  (Ref.  17) 
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to  loeata  yF  in  the  middle  of  the  first  grid  step  in  the  co-direction. 


It  would,  of  course,  be  desirable  to  obtain  analytical  expressions 
for  the  Couette-type  solutions.  Unfortunately,  however,  analytical 
solutions  of  Eqs.  (51)  to  (55)  can  be  obtained  only  after  imposing  severe 
simplifying  assumptions.  To  give  an  example,  an  analytical  solution 
of  Eq,  (51)  can  be  obtained  under  the  assumptions  (1)  that  the  exponential 
factor  in  Eq.  (56)  can  be  omitted,  (2)  that  p  =  const  (this  reduces  the 
Couette-type  solution  for  u  to  a  1 -parametric  family),  (3)  that  mw  =  0, 
(4)  that  dp/dx  =  0,  and  (5)  that  J^B  dy  =  0.  Then,  using  the  dimen¬ 
sionless  variables 


* 

y 


=  k 


y 


u 


K 


u 


(58) 


Eq.  (51)  can  be  integrated  to  obtain 


* 

u  =  - 


1  +Vl  +  4y 


?y. . -  ■,  +  In  (2y*  +  V 1  +4y*2) 


*2 


(59) 


Then,  this  solution  can  be  used  to  integrate  also  the  remaining  Eqs. 
(52)— (55)  under  the  additional  simplifying  assumptions  that  (6)  all 
laminar  and  turbulent  Prandtl  and  Schmidt  numbers  are  equal  to  1, 
and  (7)  that  the  integrals  appearing  in  Eqs.  (52)- (55)  can  be  approxi¬ 
mated  by  y  times  the  average  value  Q  of  the  integrands  in  the 
wall  region.  Using  the  dimensionless  variables 


* 

$ 


=  K 


w 


(60) 
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Eqs.  (52)  -  (55)  can  be  expressed. in  the  common  form 


=  — 
dy  !  + 


1  +  Q 


Sit 


2y 


(61) 


1  +  \i  +  4y 


*2 


which  can  be  integrated  to  yield 


+*  =  u*  ♦  °!{  VTTv^.  1 .  1„  JLlipS!?}.  (62) 


As  mentioned  above,  the  present  method  has  not  relied  on  such 
analytical  solutions,  but  has  carried  out  and  used  numerical  solutions  of 
Eqs.  (51)  -  (55).  Occasionally,  it  has  been  possible  to  fit  these  numerical 
solutions  to  suitable  analytical  formulas,  as  suggested  in  Ref.  17. 

2.  2.  5.  Matching  of  winnerw  and  "outer”  solutions 

The  first  control  volume  away  from  the  wall  utilized  for  the 
discretization  of  the  differential  equation  of  the  problem  starts  at  the 
half-way  point  co  =  ^  h^  of  the  first  grid  size  in  the  co-direction  (see 
Fig.  5).  For  the  discretization  the  dependent  variables  have  been 
assumed  to  vary  linearly  in  the  co-direction  between  grid  points.  How¬ 
ever,  since  the  first  control  volume  employed  for  this  purpose  starts 
at  co  =  -jhj,  only  the  variation  of  each  unknown  between  co  =  h^  and 
Co  =  hj  has  been  actually  used.  If  this  linear  variation  were  extrapo¬ 
lated  to  the  wall  (co  =  0)  there  could  be  found  a  value  for  each  depen¬ 
dent  variable,  called  a  "slip  value”  <t>  =  (u  ,  T  ,  c  ,  T  ),  that  could 

Ts  s  s’  a,  s  e,  s 

be  very  different  from  the  actual  value  of  the  dependent  variables 
at  the  wall  (see  Fig.  5).  This  is  due  to  the  sharp  variation  of  the  vari¬ 
ables  near  the  wall,  which  makes  the  assumption  of  linear  variation  a 
very  poor  one  in  that  region.  In  fact  care  must  be  taken  that  the  discre¬ 
tization  for  the  ’’outer”  solution  remains  valid  as  close  to  the  wall  as 
co  =  -jhj.  The  fictitious  "slip  values”  of  the  dependent  variables  are, 
however,  the  proper  boundary  conditions  to  be  used  for  the  ”outer” 
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solution.  This  means  that  the  two-dimensional  calculation  can  be  for¬ 
mally  started  "at  the  wall"  provided  the  defined  slip-values  6  (which 
will  have  to  be  found  via  the  matching  procedure)  are  used  as  boundary 
conditions  instead  of  the  physical  wall  conditions  4>w  of  the  problem. 

In  other  words,  the  presence  of  the  inner  layer  is  formally  equivalent 
to  the  introduction  of  appropriate  slip-conditions  at  the  wall. 

The  matching  procedure  by  which  the  slip  values  are  determined 

is  shown  on  Fig.  5:  The  linear  variation  of  each  dependent  variable 

between  the  slip  value  at  the  wall  and  its  value  at  the  point  w  =  h. 

1  1 

must  be  such  that  both  its  value  and  its  derivative  at  to  =  h^ 
agree  with  the  value  and  slope  of  the  Couette-type  solution  at  that  point. 
(In  a  generalization  that  has  been  found  to  improve  the  accuracy  of  the 
solution,  the  slopes  are  required  to  be  matched  at  to  =  £h^,  where 

7  *«<»■> 

Analytically,  this  matching  procedure  can  be  defined  as  follows, 
using  the  expression  given  by  Eq.  (57)  for  the  Couette-type  solution  . 
for  the  dependent  variables: 


Eqs.  (63)  constitute  a  system  of  two  equations  for  expressing  the  two 

unknowns  d>  and  F  in  terms  of  6  and  «b0.  The  first  of  these 
Ts  w  '  w  T  c- 

results,  namely 

4>s  =  (64) 

will  be  used  as  the  wall  boundary  condition  for  the  outer  solution,  while 
the  second,  namely 

F  =F  (X>,’X  )  (65) 

w  w'Y2*  Yw  '  ' 
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will  provide  the  values  of  the  wall  fluxes  after  the  outer  solution 
(and  hence  cj^)  have  been  determined. 

Alternatively,  in  case  rather  than  ^  have  been  specified, 

Eqs.  (63)  will  be  solved  for  <i>B  and  in  terms  of  F^  and 


*8  =  (62*  Fw}  (66) 

*w=  *w(*2’  Fw}  (67) 


Finally,  note  that  the  boundary  conditions  expressed  by 
Eq.  (64),  or  Eq.  (66),  refer  to  the  unknowns  cj>,  which  do  not  include 
r.  A  boundary  condition  for  the  latter  is  also  needed  in  order  to  solve 
the  system  of  equations  (49).  This  conditions  is  obtained  simply  by 


matching  the  linear  variation  of  t  between  t 

(namely 


and  with  the  value 


of  the  Couette-type  variation  of  r 
1  , 

"■  r  V 


T  =  fl 


3u 
T  c)y 


)  at  the  point 


t2+ts 


(68) 


du 


Note  that  the  Couette-flow  solution  for  u,  and  hence  for  t  =  ♦  has 

already  been  obtained  as  the  first  component  of  Eq.  (64),  or  Eq.  (66). 


2.  2.  5.  Solution  by  elimination 

The  system  of  finite -difference  equations  (4.9),  which  can  be 
written  in  the  form 


T 


+  B.T.  , 

i  l-l 


(i  =  2,  3,  .  .  . ,  n-1) 


(69) 


together  with  the  boundary  conditions  constitute  a  system  of  n 
equations  for  the  n  unknowns  =  (u.,  t.,  T.,  c  .,  T  .),  where  the 

X  XXX  Qj  X  0 1  X 

subscript  i  corresponds  to  each  grid  point  in  the  co-direction.  Note 
that  each  equation  (69)  is  the  result  of  discretization  in  the  control 
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volume  centered  around  the  point  i.  As  pointed  out  earlier,  the 
first  control  volume  away  from  the  wall  that  is  used  for  discretization 
of  the  differential  equations  is  the  one  centered  around  the  point  2, 
while  the  last  is  the  one  centered  around  the  point  n-1.  This  is 
reflected  by  the  range  of  i  in  Eq.  (69).  The  two  boundary  conditions 
have  the  same  general  form  as  Eqs.  (69)  and  can  therefore  be  con¬ 
sidered  as  the  finite -difference  equations  corresponding  to  points  i  =  1 
and  i  =  n.  As  already  mentioned,  the  one  corresponding  to  i  =  1  is 
made  up  of  Eqs.  (64)  -  or  (66)  -  together  with  Eq.  (68)  and  gives  the 

values  of  7  =  in  terms  of  the  values  of  IT.  and  ^  (or  IT  ): 

si  2  w  w 

the  one  corresponding  to  i  =  n  simply  fixes  the  values  of  ^  by 
equating  them  to  the  free -stream  values. 

Since  the  value  of  "3T  (or  F  )  is  given,  the  inner  boundary 
condition  (i  =  i)  can  be  written  in  the  form 

^l=V^2  +  Ci  (70) 


Then,  by  substituting  in  Eq.  (69)  for  i  -  2,  we  find 

VA2*  V^2  <71> 

with  A'2  =  (1  -  B2-  A1)'1A2  and  =  (1  -  Ey  Al)"l(B2*  Cfj  +  cTj).  Con¬ 
tinuing  this  process,  all  3-term  recursion  equations  (69)  can  be  reduced 
to  2-point  recursion  equations  of  the  form 


i 


+  c. 

l+i  i 


(72) 


(i  =  2,  3, ....  n-1) 


where  the  matrix  a|  and  the  vector  c!  are  found  by  substituting  ^ 
as  given  by  Eq.  (72)  in  Eq.  (69): 


l 


+  B.  •  (A.  •  ®.  +  C  .  .)  +  C. 

i  '  l-l  l  1-1'  i 


This  results,  for  i  =  2,  3,  .  .  . ,  n-1,  in  the  recursion  relations 
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a!  =  (i  -  b.  •  a!  ,)_1  a. 

i  '  i  l-i7  i 

C.  =  (i  -  B.-  a!  j)-1.  (B.  -  c!  ^C.) 


It  is  clear  from  Eq.  (70)  that  the  coefficients  Aj  =  A^  and 
C  j  =  Cj  are  given  by  the  inner  boundary  condition.  Consequently, 
starting  with  i  =  2,  Eqs.  (73)  can  be  applied  for  each  subsequent  i 
to  yield  to  actual  numerical  values  for  all  coefficients  a!  and  C  ! 
up  to  i  =  n- 1. 


When  A^  ^  and  ^  are  known,  then  Eq.  (72)  for  i  =  n-1, 

namely  the  relation 


=  A 


n' 


.  9  +  C  . 

■  In  n-1 


(74) 


yields  the  numerical  value  of  sTn  since  the  value  of  is  given  as 
the  outer  boundary  condition. 


Then,  using  Eqs.  (72)  successively  in  the  order  i  =  n-2,  n-3, 

. . . ,  2  together  with  the  computed  coefficients  a!  and  (T.,  the  numeri¬ 
cal  values  of  all  unknowns  ®.  inside  the  boundary  layer  are  found  at  the 
station  x  =  X,.  Note  that  the  slip  value  IT  is  also  available  through 
Eq.  (70). 

The  elimination  procedure  described  above  is  analogous  to 
the  more  elaborate  direct  method  used  for  the  solution  of  the 
streamf unction  equation  in  the  electrical  problem  (see  Ref.  1).  In 
fact,  the  method  described  above  is  equivalent  to  considering  the  un¬ 
known  $  .  asa  "fundamental  unknown"  in  the  context  of  the  termin- 

n-i 

ology  used  in  Ref.  1. 
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2.  2.  6,  Definition  of  the  external  boundary  and  of  the  entrainment 
rate 

When  considering  flows  where  the  external  boundary  is  a  free 
boundary,  there  is  considerable  freedom  in  defining  the  position  yE(x) 
of  this  boundary  for  mathematical  convenience.  There  is  of  course  a 
one-to-one  correspondence  between  the  criterion  selected  for  the 
definition  of  the  external  boundary  and  the  resulting  entrainment  rate 
mg  through  this  boundary.  Note  that  the  entrainment  rate  enters  in 
the  solution  of  the  boundary  layer  problem  both  by  appearing  explicitly 
in  the  differential  equations  through  the  quantity  R  [see  Eq.  (40a)] 
and  by  defining  the  normalization  of  the  dimensionless  cross-stream 
variable  u. 

An  advantageous  choice  of  the  E-boundary  is  to  define  it  as  the 

line  on  which  the  streamwise  velocity  u  equals  a  given  percentage  of 

the  free.- stream  value  U  .  This  choice  has  already  been  introduced 

oo 

in  the  process  of  formulating  this  model  (see  Section  II).  The  definition 
yE(x)  =  6ggg(x)  is  consistent  with  the  definition  used  to  obtain  measure¬ 
ments  of  the  entrainment  rate  (Refs.  16  and  22).  These  measurements 
were  used  to  derive  numerical  values  of  the  turbulent  structure  param¬ 
eter  a^  through  Eq.  (26). 

2.  2.  7.  Turbulence  structure  parameters  a^f  a^,  a^ 

Values  for  the  turbulence  structure  parameters  a^,  a ^ 
measured  and  reported  by  Bradshaw  (Ref.  16  and  22)  have  been  used 
in  carrying  out  the  numerical  computations  reported  herein. 

These  measurements  were  carried  out  under  a  variety  of  boundary 
layer  conditions  and  constitute  limited  experimental  proof  of  the  uni¬ 
versality  of  these  parameters.  The  parameter  a^  has  thus  been 
shown  to  be  constant  across  the  boundary  layer,  with  a  numerical 
value  of  0.  15;  the  numerical  values  of  a^  and  a^  used  in  the  compu¬ 
tations  are  given  in  Table,  1  as  functions  of  the  dimensionless  boundary 
layer  coordinate  y/6. 
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TABLE  1 


a2  X  102 


a 


3 


0.  0 

0.  0 

0.  0 

0.  1 

4.  0 

0.  2 

0.2 

7.  5 

1.0 

0.  3 

8.  8 

2.  0 

0.4 

9.  2 

3.  0 

0.  5 

9.  5 

5.  0 

0.  6 

9.  2 

7.  7 

0.  7 

8.  5 

14.  0 

0.  8 

7.  4 

23.  5 

0.  9 

6.  0 

30.  0 

1.  0 

4.  1 

33.  0 

1.  1 

2.  0 

36.  0 

1.  2 

0.  5 

37.  3 

1.  3 

0.  1 

38.  5 

1.4 

0.  0 

41.  5 
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Appendix  C  describes  the  computer  program  LAYER,  which  has 
been  coded  in  FORTRAN  for  the  solution  of  the  gasdynamic  problem  in  the 
compressible,  turbulent  magnetohydrodynamic  boundary  layer.  It  includes 
a  list  of  routines  and  important  variables,  specification  of  required  input 
data  and  formats,  and  a  description  of  available  output. 
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IV.  APPLICATIONS 

As  the  solutions  to  each  of  the  subproblems  described  in  Section  1 
were  obtained,  appropriate  test  cases  were  analyzed,  with  the  dual  intent 
of  checking  the  solution  procedure  and  developing  insight  into  the  physical 
roles  of  the  mechanisms  under  study. 

Subsequent  to  the  separate  development  work  on  each  subprob¬ 
lem,  the  electrical  solution  over  the  entrance  region  and  the  solution  of 
the  gasdynamic  problem  were  coupled  together  and  used  to  analyze  the 
electrical  characteristics  of  the  Hirho  accelerator  at  AEDC. 

In  this  Section  the  results  obtained  by  solving  each  subproblem 
in  various  geometries  and  operating  conditions  will  be  described. 


1.  SOLUTIONS  OF  THE  ELECTRICAL  PROBLEM 

To  apply  the  solution  of  the  electrical  problem  over  the  entire 
channel,  it  was  natural  to  select  cases  that  had  been  investigated  exten¬ 
sively  during  the  previous  two  phases  of  this  study.  Accordingly,  the  first 
test  case  chosen  was  that  of  a  J  X  B  accelerator  as  analyzed  in  Refs.  1 

and  2  operating  with  nitrogen  seeded  with  0.  5%  (by  weight)  potassium  at 

2 

p  =  0.  5  atm,  with  average  J  =30  A/cm  and  uniform  magnetic  induction 

/  2  y 

of  2  Wb/m  .  The  geometry  of  each  electrode  pair  was  as  in  Refs.  1  and 

2,  namely  12  mm  for  the  conductor  and  3  mm  for  the  insulator  segment, 
with  an  interelectrode  distance  of  20  mm.  For  the  entrance  region 

it  was  found  sufficient  to  consider  six  electrode  pairs  (of  identical 

dimensions  i  .  =  12  mm  and  f  j  .  =  3  mm)  preceded  by  57  mm  of 
c,  l  ^  i 

insulating  wall.  As  in  Ref.  2,  three  computations  were  carried  out, 
differing  only  in  the  assumed  T  and  U  profiles:  the  first  computation 
used  uniform  gas  temperature  T  =  3000  °K  and  uniform  gas  velocity 
U  =  3000  m/s;  the  second  computation  used  uniform  T  =  3000  °K  and 
uniform  gas  velocity  300  m/s,  one  order  of  magnitude  smaller  than  in 
the  first  computation;  the  third  computation  accounted  for  the  variation 
of  both  T  and  U  in  the  electrode  wall  boundary  layers,  as  described 
in  Ref.  2. 
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Fig,  6  presents  the  current  distribution  and  Fig.  7  the  electric 
potential  distribution  corresponding  to  each  of  these  computations. 

These  computations  illustrate  again,  as  did  the  computations 
described  in  Ref.  2,  how  the  local  electrical  behavior  is  influenced  by 
finite  reaction  rates,  electron  energy  convection  and  thermal  and 
velocity  boundary  layers.  These  effects  have  been  discussed  in  Ref.  2. 

In  addition,  however,  the  present  solutions  over  the  entire  entrance 
region,  illustrate  the  following  facts: 

(1)  The  physical  extent  of  current  fringing  upstream  of  the 
first  electrode  is  small  in  this  seeded  nitrogen  accelerator,  not  exceeding 
approximately  one  electrode  period  L.  This  shows  that  the  finite  insu¬ 
lator  length  Lq  used  upstream  of  the  first  electrode  is  more  than  suffi¬ 
cient  for  the  mathematical  treatment  of  the  entrance  region  of  this 
accelerator. 

(2)  The  current  and  electric  field  distributions  (see  Figs.  6  and 
7),  as  well  as  the  distributions  of  n  ,  T  and  of  all  plasma  properties, 
can  be  considered  periodic,  to  a  very  good  approximation,  after  the  third 
or  fourth  electrode  pair.  This  shows, again,  that  the  number  of  electrode 
pairs  used  for  the  mathematical  treatment  of  the  entrance  region  is 
more  than  sufficient  for  this  accelerator.  Comparison  of  the  periodic 
parts  of  Fig.  6  with  the  corresponding  current  distributions  computed  in 
Phase  II  under  the  imposed  assumption  of  periodicity  (Figs.  1,  4,  7  of 
Ref.  2)  show  detailed  agreement. 

The  limited  current  fringing  and  rapid  approach  to  periodicity 
exhibited  by  these  computations  is,  of  course,  dependent  upon  the  choice 
of  the  working  fluid.  On  the  other  hand,  computations  performed  by  this 
method  for  monatomic  seeded  and  unseeded  gases  have  indicated  similar 
characteristics:  current  fringing  extended  at  most  a  few  electrode 
periods  upstream  of  the  first  electrode  pair,  while  periodicity  was 
reached  by,  at  most,  the  eighth  electrode  pair. 

Finally ,  the se  computations  have  illustrate'’  that  current  fringing 
in  the  exit  region  of  J  X  B  accelerators  is  not  more  severe  than  the 
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Fig.  6.  Current  distributions  in  seeded-nitrogen  accelerator 
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Fig.  7.  Electric  potential  distributions  in  seeded-nitrogen  accelerator 


Uniform  T  =  3000  °K  and  U  =  300  m/s 


Boundary  layer  variation  of  T  and  U  included 


AEDC-TR  70-86 


AEOC-TR-70-86 


fringing  reported  above  for  the  entrance  region. 

2.  BOUNDARY  LAYER  SOLUTIONS 

The  novel  formulation  and  method  of  solution  used  in  this  study 
for  the  boundary  layer  problem  has  provided  inducement  for  careful 
testing  of  their  accuracy  and  reliability.  For  this  purpose,  consistent  and 
well  documented  cases  of  turbulent  boundary  layer  flows  were  sought, 
for  which  experimental  results  were  also  available.  It  was  natural, 
therefore,  to  select  test  cases  from  the  33  incompressible  turbulent 
flows  documented  by  D.  Coles  for  the  AFOSR-IFP-Stanford  1968 
Conference  on  Turbulent  Boundary  Layers  (Refs.  23  and  24).  The 
experimental  results  collected  for  this  conference  (Ref.  24)  and  the 
results  of  the  computations  carried  out  by  the  participants  (Ref.  23) 
concerned  the  axial  variation  of  the  skin  friction  coefficient  c„  the 
shape  factor  H  =  6  /G,  and  the  momentum-thickness -Reynolds  number 
Re^=  U^G/v,  where  5  is  the  displacement  thickness,  G  the  momen¬ 
tum  thickness  and  v  the  kinematic  viscosity.  These  particular  quanti¬ 
ties  were  selected  on  the  basis  that  they  should  provide  a  more  sensitive 
indication  of  the  accuracy  of  the  prediction  method  than  the  overall 
velocity  profile.  For  the  purpose  of  the  conference,  values  of  these 
quantities  were  obtained  numerically  from  measured  velocity  profiles 
at  various  axial  stations.  Specifically,  c^  was  obtained  by  fitting  the 
"logarithmic  law  of  the  wall"  to  the  experimental  velocity  profile,  and 

6  and  G  were  calculated  using  a  modified  Simpson  integration.  The 

1 

predictions  reported  in  Ref.  23  were  all  started  at  a  specified  station  in 
each  case,  with  initial  conditions  such  as  to  match  c^  and  Re^  with  the 
experimental  values  at  that  station. 

Two  of  the  cases  that  were  selected  for  testing  of  the  present 
method,  and  which  will  be  described  here,  were  those  identified  in 
the  conference  as  cases  2100  and  2400.  Case  2100  was  subject  to  an 
initial  mildly  negative  pressure  gradient,  which  then  became  strongly 
positive  at  about  x  =  18  feet,  leading  to  eventual  separation.  Case  2400 
involved  a  moderately  positive  pressure  gradient,  which  dropped  abruptly 
to  zero  near  x  =  5  feet. 
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The  numerical  results  obtained  for  these  two  cases  by  the  various 
participants  in  the  conference  are  shown  on  Figs.  8  and  12  respectively. 
(These  figures  are  reproduced  here  from  Ref.  23.  )  The  results  of  the 
new  solution  method  developed  in  this  study  are  shown  on  Figs.  9-11  for 
Case  2100,  and  on  Figs.  13-15  for  Case  2400;  note  that  the  latter  results 
were  obtained  by  matching  the  initial  values  of  C£  and  of  the  kinematic 
viscosity  only. 

It  is  clear  from  Figs.  8-15  that  the  new  method  of  solution 
developed  in  this  study  gives  quite  satisfactory  results  for  these  two 
incompressible  test  cases,  and  in  fact  shows  up  quite  well  in  comparison 
to  conventional  methods.  In  particular,  the  results  for  C£  are  in  very 
good  agreement  with  the  experimental  data. 

Results  obtained  by  the  new  method  of  solution  for  the  compres¬ 
sible  magnetohydrodynamic  boundary  layer  development  in  the  Hirho 
experiment  will  be  described  in  the  following  subsection. 

3.  COMPUTATIONS  PERTAINING  TO  THE  HIRHO  EXPERIMENTS 

The  complete  two-dimensional  solution  of  the  coupled  electrical 
and  gasdynamic  problems  has  been  carried  out  over  the  entire  powered 
section  of  the  AEDC  Hirho  channel  (Ref.  8)  with  seeded  air  as  the 
operating  fluid.  This  was  the  first  application  of  the  complete  method; 
it  included  two-dimensional  solution  of  the  electrical  problem,  quasi- 
one -dimensional  solution  of  the  gasdynamic  problem  in  the  core  of  the 
flow,  and  two-dimensional  boundary -layer  solution  on  the  electrode 
walls.  This  application  gave  impressive  demonstration  of  how  essential 
it  is  to  avoid  oversimplifications  when  attempting  to  predict  or  inter¬ 
pret  the  behavior  of  multielectrode  J  X  B  devices.  It  has  already  been 
found  for  the  Hirho  channel  (Ref.  8),  as  for  many  other  cases,  that 
as  long  as  one -dimensional  analyses  are  employed  there  results  con¬ 
siderable  discrepancy  between  predicted  and  observed  Hall  fields.  This 
work  has  shown  that  correct  prediction  and  interpretation  of  the  Hall 
field,  which  is  an  important  design  parameter  for  magnetohydrodynamic 
channels*  cannot  be  attained  even  by  the  two-dimensional  but  linearized, 
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Fig.  13.  vs.  x  for  the  conditions  of  Case  2400. 
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"constant-property,  "  computation. 

The  present  method  has  been  very  successful  in  predicting  and 
interpreting  the  experimental  results  (Ref.  8)  that  have  been  obtained 
in  the  Hirho  channel.  In  particular,  it  demonstrated  the  substantial 
effect  of  the  local  axial  current  density  on  the  Hall  field,  especially 
at  relatively  low  values  of  the  Hall  parameter  (3/e.  The  presence  of 
finite  (and  negative)  in  the  core  of  J  X  B  accelerators  tends  to  re¬ 
duce  the  Hall  field  Ex,  even  when  there  is  no  net  axial  current  leakage, 
as  can  be  clearly  seen  by  writing  Ohm's  law  in  the  x-direction  in  the  form 

E  +  K  =  —  J  [  -  tg  p  +  |  ]  (75) 

where  K  the  thermal  diffusion  vector  (Refs.  1  and  2)  and  <p  the  angle 
the  current  density  J  makes  with  the  positive  y-axis.  First,  a  finite 
Jx  leads  to  finite  tg  <p  which  diminishes  the  absolute  value  of  the  square 
bracket  in  Eq.  (75).  Secondly,  and  more  importantly,  a  finite  Jx 
increases  the  rate  of  local  Ohmic  heating  and  can  thus  result  in  sub¬ 
stantially  higher  values  of  the  scalar  conductivity  <r,  which  appears  in 
the  denominator  in  Eq.  (75). 

The  cumulative  effects  of  these  two  mechanisms  on  the  Hall  field 
in  the  Hirho  channel,  as  investigated  by  this  study,  are  illustrated  on  Fig. 
16:  The  quantity  chosen  for  this  illustration  is  the  potential  of  each 
cathode  in  the  channel  relative  to  cathode  No.  1;  this  choice  allows 
direct  comparison  with  experimental  measurements  (shown  as  curve  1 
on  Fig.  16).  The  operating  conditions  were  those  of  Hirho  Run  No. 

1412  (Ref.  8).  Curve  4  shows  the  results  of  a  one -dimensional  compu¬ 
tation,  performed  by  the  authors  of  Ref.  8,  which  shall  be  called  compu¬ 
tation  I.  Computation  I  estimated  the  average  J  from  the  given  elec¬ 
trode  currents,  and  assumed  Jx  to  be  zero  (or  at  most  a  small  fraction 
of  Jy,  estimated  as  in  Ref.  25,  but  without  any  effect  on  the  local  Ohmic 
heating).  Curve  3  shows  the  results  of  a  two-dimensional  solution  of  the 
electrical  problem  by  the  present  method,  which  however  used  as  inputs 
the  distributions  of  the  gasdynamic  unknowns  U  and  T  obtained  from 
computation  I.  Finally,  curve  2  was  obtained  by  coupling  the  electrical 
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volts 


Fig.  16.  Cathode  potentials  relative  to  cathode  #1  for  a  J  X  B 
accelerator,  [  Operating:  conditions  correspond  to 
Hirho  run  1412  {Ref.  8).  ] 
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and  gasdynamic  problems  over  the  first, four  electrodes,  and  accounting 
fully  for  the  increased  Ohmic  heating  due  to  finite  local  J.  Curve  2 
demonstrates  the  ability  of  this  method  to  predict  and  interpret  the 
operating  characteristics  of  multielectrode  J  X  B  devices. 

After  the  effect  of  finite  J^was  demonstrated,  additional  compu¬ 
tations  were  performed  to  investigate  possibilities  of  improving  the 
performance  by  diminishing  the  angle  <p  through  changes  in  the 
electrode  geometry.  Thus  a  second  run  treated  a  "staggered-electrode" 
geometry,  where  the  cathodes  were  shifted  in  the  downstream  direction 
by  5/9  of  the  electrode  period  (i.  e. ,  approximately  l")  compared  to 
the  position  of  the  corresponding  anodes,  and  a  third  run  treated  a 
geometry  with  finer  electrode  segmentation  (L,/ D  =  5/9  instead  of  9/8). 
The  operating  conditions  used  in  all  computations  were  those  of  Run  No. 
1412  (Ref.  8).  The  results  are  shown  on  FigB.  17  and  18. 

Figure  17(a)  shows  the  current  distribution  in  the  actual  Hirho 
geometry  as  computed  by  the  idealized  "constant  property"  method,  to 
contrast  with  Fig.  17(b)  which  shows  the  current  distribution  as  computed 
by  the  present  method,  including  thermal  diffusion,  finite  reaction  rates, 
and  electron  energy  convection.  The  generally  greater  angle  <p  between 
the  current  lines  and  the  y-axis  (and  correspondingly  greater  Jx  for 
given  Jy)  is  the  most  noticeable  characteristic  of  Fig.  17(b)  as  com¬ 
pared  to  Fig.  17(a).  The  potential  distribution  corresponding  to  the 
current  distribution  of  Fig.  17(b)  is  shown  on  Fig.  17(c),  on  which  the 
contours  are  equipotentials  for  the  electric  field  E  (unprimed)  and 
the  contour  interval  is  50  volts. 

Figure  18(a)  shows  the  current  distribution  in  the  "staggered- 
electrode"  geometry,  and  demonstrates  that  in  this  case  the  angle  <p 
has  an  average  value  of  the  order  of  5°,  as  opposed  to  20°  in  the  actual 
geometry  (Fig.  17(b)).  Finally,  Fig.  18(b)  shows  the  current  distri¬ 
bution  in  the  case  of  finer  segmentation  and  demonstrates  that  in  this 
case  the  average  value  of  the  angle  q>  is  of  the  order  of  10°. 
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Fig.  17.  Current  and  electric  potential  distributions  for  actual 
Hirho  geometry. 


(a)  Current  distribution  according  to  "constant  property"  computation 


(b)  Current  distribution  according  to  finite  rate  computation 


(c)  Electric  potential  distribution  according  to  finite  rate  computation 
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Fig.  18.  Current  distributions  corresponding  to 
modified  Hirho  geometries. 


(a)  Current  distribution  for  staggered  geometry 


(b)  Current  distribution  for  more  finely  segmented  electrodes 
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Finally,  the  application  of  the  two-dimensional  boundary -layer 
solution  gave  an  illustration  of  the  boundary-layer  development  in  the 
Hirho  channel  under  the  operating  conditions  of  Run  No.  1412  (Ref.  81. 
The  reader  is  reminded  that  the  working  fluid  was  air  seeded  with  0.  2% 
(by  weight)  potassium  and  that  the  flow  conditions  at  the  selected  initial 
station  x  =  0  (at  d  =  44.  88  cm  from  the  nozzle  throat  —  see  Ref.  8) 
were  as  follows: 

U  =  2947  m/s 
co  ' 

T  =  3017  °K 
00 

T  =  300  °K 
w 

p  =  4.  56  atm 
=  -  2  X  106  N/m3 
Re  =  U  d/v  =  8.  7  X  106 
6  =  6.  1  mm 

<y  =  4.  07  X  105  A/m2 
B  =  2.  85  Wb/m2 

For  the  initial  profiles  at  x  =  0  a  i/7th  power  law  variation  has  been 

assumed  for  the  mean  velocity  u  and  for  the  temperature  ratio 

(T-T  )/(T  -T  ).  The  axial  variations  of  U  and  of  the  current 
'  w  '  oo  w  oo 

density  components  Jx  and  J^.  were  obtained  from  the  quasi -one - 
dimensional  solution  of  the  core  problem  and  the  two-dimensional  solu¬ 
tion  of  the  electrical  problem,  and  have  been  expressed  as  follows: 


U  (x)  =  2947  +  2250x 
oo 

m/  s 

Jx(x)  =  -  (2.  52  +  5.  96x)  ■  105 

A/m 

Jy(x)  =  (4.  07  +  34.  6x)  ■  105 

A/m 

The  solution  has  provided,  in  addition  to  the  description  of  boundary- 
layer  development,  the  distributions  of  skin  friction  and  wall  heat 
transfer  rate  along  the  channel.  The  results  are  shown  on 
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Figs.  19-25.  Figures  19,  20,  21,  22,  and  23  present  the  profiles  of  the 
unknowns  u,  t,  T,  n  ,  and  T  respectively,  at  various  x-stations. 

Note  that  Fig.  20  shows  the  profile  of  the  turbulent  shear  stress  t 
alone.  If  the  laminar  contribution  were  added,  the  total  shear  stress 
Tj^+t  would  have  a  maximum  at  the  wall  due  to  the  effective  favorable 
pressure  gradient  in  the  supersonic  diverging  channel.  Fig.  21  shows 
how  the  mean  gas  temperature  profile,  which  was  initially  assumed  to 
follow  a  l/7  power  law,  is  affected  by  Ohmic  heating  and  dissipation  so 
as  to  develop  a  peak  close  to  the  wall  at  downstream  stations.  It  should 
be  pointed  out  that  a  more  realistic  initial  profile  could  be  prescribed 
in  terms  of  the  turbulent  version  of  Crocco's  relation  (Ref.  26) 


c  T  +  iu2 
P  2 


*  ^[(1-^ 


u 

XT  ' 

00 


T,] 


(76) 


where  is  the  freestream  stagnation  temperature  and 


T  /T°  . 
w'  00 


The  computed  results  for  the  variation  of  the  skin  friction  coef- 

0T 

ficient  c,  and  the  heat  transfer  coefficient  c„  =  -(k  £i)  /[  p  U  (H  -H  )] 

f  H  oy  w  r  00  00  00  w 

along  the  channel  are  shown  on  Fig.  24.  The  skin  friction  results  lie 

between  those  obtained  by  the  method  of  Enkenhus  and  Mahez  and  by  that 

of  Elliot,  Bartz,  and  Silver  and  reported  in  Ref.  8.  The  computed  growth 

of  the  boundary  layer  thickness  6  and  of  the  displacement  thickness  6^ 

are  shown  on  Fig.  25.  The  latter,  by  definition,  reduces  the  effective 

channel  height  D  for  the  core  flow.  Taking  into  consideration  that  the 

channel  height  D  at  the  initial  station  x  =  0  is  3.  28  cm,  and  the  angle 

of  divergence  of  the  channel  is  1  26  ,  we  see  that  the  ratio  6  /D  varies 

-3  -2 

from  approximately  6XiO  at  x  =  0  to  3  X  10  at  x  =  14  cm. 
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Fig.  19.  Development  of  mean  velocity  profile  in  a  JXB  accelerator 
[Operating  conditions  correspond  to  Hirho  run  1412  (Ref.  8).  ] 


Fig.  20.  Development  of  turbulent  Bhear  stress  profile  in  a  J  X  B  accelerator 
[Operating  conditions  correspond  to  HirHo  run  1412  (Ref.  8).  ] 
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Fig.  22.  Development  of  mean  electron  number  density  profile  in  a  J  X  B  accelerator 
[Operating  conditions  correspond  to  Hirho  run  1412  (Ref.  8).  ] 
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Fig.  23.  Development  of  mean  electron  temperature  profile  in  a  J  X  B  accelerator 
['Operating  conditions  correspond  to  Hirho  run  1412  (Ref.  8),  ] 
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Fig.  24.  Streamwise  variation  of  skin  friction  and  wall  heat  transfer 
rate  coefficients  in  a  J  X  B  accelerator.  [  Operating 
conditions  correspond  to  Hirho  run  1412  (Ref.  8).  ] 
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Fig.  25.  Growth  of  boundary  layer  ^hickness  6  =  6995  and 
displacement  thickness  6  .  [  Operating  conditions 

correspond  to  Hirho  run  1412  (Ref.  8),  J 
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APPENDIX  A 

DESCRIPTION  OF  COMPUTER  PROGRAM  "INLET" 

A.  i.  GENERAL  DESCRIPTION 

Program  INLET  has  been  written  to  solve  the  coupled  electrical 
+  gasdynamic  problem  over  the  whole  inlet  or  exit  regions  (see  main 
text).  The  user  may,  at  his  option,  specify  solution  of  the  electrical 
part  of  the  problem  alone,  or  of  the  coupled  problem;  in  the  latter  case, 
the  gas  velocity  and  temperature  variations  are  calculated  along  the 
streamlines,  without  viscosity  or  heat  conduction  effects.  The  solution 
may  be  performed  over  any  specified  geometry  incorporating  up  to 
twenty  electrode  pairs  (powered  and  unpowered)  bounded  by  insulating 
walls  or  insulator  segments  at  the  upstream  and  downstream  ends. 

INLET  can  treat  any  working  fluid,  distinguishing  as  many  as  25  sepa¬ 
rate  components,  formed  from  up  to  10  distinct  elements  and  taking 
part  in  up  to  22  reactions.  Channel  geometry  (including  dimensions  of 
all  conductor  and  insulator  segments),  axial  distribution  of  magnetic 
field  strength,  electric  current  passed  through  each  electrode  pair, 
distribution  of  streamfunction  across  the  channel  at  upstream  and  down¬ 
stream  ends  (when  required),  and  initial  profiles  of  gasdynamic  un¬ 
knowns  must  all  be  provided.  INLET  will  then  proceed  for  a  specified 
number  of  iterations  between  gasdynamic  and  electrical  solutions, 
whereupon  a  restart  deck  may  be  generated  in  order  to  continue  the 
calculation  from  this  point  at  a  later  time. 

INLET  has  been  coded  in  FORTRAN  for  Control  Data  Corporation 
6000  series  computers,  requiring  approximately  100,  000  memory  loca¬ 
tions  and  capacity  to  retain  28  significant  digits  (double  precision). 
Typical  calculations  require  approximately  40  seconds  of  central  pro¬ 
cessor  time  on  the  CDC  6600  for  each  complete  (electrical-gasdynamic) 
iteration  cycle. 

The  following  sections  provide  a  more  detailed  description  of 
program  INLET.  Section  A.  2  describes  input  data  and  formats,  A.  3 
lists  available  output,  A.  4  lists  the  routines  included  in  the  program, 
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and  Section  A.  5  defines  important  variable  names  used  in  these  routines. 

A.  2.  INPUT  DATA  AND  FORMATS 
Format 

10A8  The  first  card  contains  the  HEADING  to 

be  printed  at  the  top  of  each  new  page. 
Leaving  columns  1-8  blank  terminates 
the  program. 

813  The  second  card  controls  the  OPTIONAL 

METHODS  OF  SOLUTION.  Each  field 
on  this  card  controls  one  of  the  options 
listed  below,  and  a  non-zero  punch  in 
any  field  selects  the  corresponding 


option. 

OPTION  INDEX 

Compute  electric  field  and 

potential  1 

Solve  for  gasdynamic  variables 

U,  T  2 

Use  of  implicit  algebraic  equation 

for  T  3 

e 

Include  effects  of  thermal  and 

concentration  diffusion  4 

Impose  periodicity  over  one 


electrode  period  as  boundary 
condition  on  5 

If  field  contains  1,  at  upstream 
end 

If  2,  at  downstream  end 
Assume  instantaneous  Saha 

equilibrium  6 

Assume  instantaneous  electron 

energy  relaxation  7 

Assume  (6)  and  (7)  during  first 

cycle  8 


Columns 

1-80 


1-24 
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Columns  Format 

1-3  13  Number  of  current  density  iteration  cycles 

to  be  performed  (up  to  25) 

1-30  1013  The  fourth  card  is  the  INPUT  CONTROL 

card.  The  remaining  input  is  separated 
into  nine  categories,  corresponding  to 
the  first  nine  fields  on  this  card.  To 
select  input  of  these  categories  a  non¬ 
zero  punch  is  placed  in  the  appropriate 
fields.  A  non-zero  punch  in  the  tenth 
field  indicates  that  a  computation  in 
process  has  been  interrupted  for  input 


and  is  to  be  continued. 

INPUT  CATEGORY  INDEX 

Output  control  deck  1 

Specification  of  working  fluid  2 

Specification  of  reactions  3 

Specification  of  geometry  and 

grid  4 

Specification  of  problem 

parameters  5 


Specification  of  initial  velocity, 
temperature  and  number 


density  profiles  6 

Specification  of  axial  distribution 
of  velocity  and  temperature  7 

Specification  of  initial  electron 

temperature  distribution  8 

Restart  deck  9 
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Category  1  OUTPUT  CONTROL 


The  first  8  cards  in  this  category  control  printout  of  the  current 
value s  of  quantities  of  interest  at  the  end  of  each  iteration  cycle.  Each 
card  corresponds  to  one  of  the  following  output  categories,  and  to  select 
output  of  any  category  at  the  ends  of  any  given  cycles,  punch  the  numbers 
of  those  cycles  (1  through  25)  in  successive  fields  on  the  appropriate 
card.  A  non-zero,  negative  number  punched  in  the  first  field  on  any 
card  will  select  that  output  at  the  end  of  every  cycle. 


Columns 


F  ormat 


OUTPUT  CATEGORY  INDEX 

Stream  function  field  1 

Current  density  field  2 

Electron  number  density  field  3 

Plasma  property  fields  ~  ^  4 

T  field  5 

e 

increment  (from  previous 
cycle)  field  6 

T,  U  fields  7 

Electric  field  and  potential  8 


1-75 


2513  Format  of  each  OUTPUT  CONTROL 


card 

The  next  2  cards  in  this  category  control  punched  output  (des¬ 
cribed  in  Section  B.  3),  suitable  for  computer -generated  contour  plots 
of  electron  temperature  or  number  density,  current  streamfunction,  or 
electric  potential,  or  for  continuing  the  computation  at  a  later  date. 


PUNCH  CATEGORY  INDEX 

Plot  deck  9 

Restart  deck  10 

1-75  2513  Format  of  PUNCH  CONTROL  cards 

The  last  card  in  this  category  indicates  plot  decks  desired. 

1-12  413  Quantities  for  which  plot  decks  are 

desired  are  selected  by  a  non-zero 
punch  in  the  appropriate  field  on 
this  card 
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Quantity 


Format 


QUANTITY 
Electron  temperature 
Electron  number  density 
Current  streamfunction 
Electric  potential 


Category  2  WORKING  FLUID 


Ratio  of  specific  heats 


(c  /c 

v  p/  V 

Number  of  elements  (up  to  10) 


1-10  E10.3 

1-5  15 

6-10  15  Number  of  components  (up  to  25) 

Each  element  is  described  on  a  card  of  the  following  format 
1-10  A10  Name  of  element 

11-20  E10.  3  Atomic  weight  of  element  (amu) 

Each  component  is  then  described  on  a  card  of  the  following 


format 

1-4 

11-60 


A4  Name  of  component 

10A5  The  composition  of  this  component  in 

terms  of  the  specified  elements  is 
defined  by  punching,  in  the  appropriate 
fields,  the  number  of  particles  of  the 
corresponding  elements  contained  in 
each  particle  of  this  component. 

1-10  E10.  3  Excitation  energy  of  each  component 

above  its  neutral  ground  state  (eV) 

Integrated  electron-neutral  collision  cross  sections,  energy  loss 
factors,  and  weighting  factors  are  given  for  each  neutral  component  in 
tabular  form  as  functions  of  electron  temperature. 

1-10  E10.3 

11-20  E 10.  3 

21-30  E10.  3 

Each  table  is  input  in  the  following  format,  where  more  than  one 
card  will  generally  be  required. 


Minimum  value  of  T  for  tables 

e 

Tabular  interval  in  T  for  tables 

e 

Maximum  value  of  T  for  tables 

e 


(°K) 

(°K) 

(°K) 
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Columns 

Format 

1-18 

18x 

19-28 

E10.  3 

29-38 

E10.  3 

39-48 

E10.  3 

49-58 

E10.  3 

May  be  used  for  identification 
Electron-neutral  collision  cross  section 

Q  (m2) 

Weighting  factor 
Weighting  factor  A^ 

Electron -neutral  energy-loss  factor  6^ 


Category  3  CHEMICAL  REACTIONS 


1-5 


15 


Number  of  reactions  (up  to  22) 


Each  reaction,  which  may  be  written  symbolically  as 


4 

Z  v-a- 

pi  J  J 


k' 

r 


4 

V 

L 

3  =  1 


i  t 

v .  a.. 
J  J 


is  characterized  by  specifying  the  components  a.,  a.  taking  part,  their 

i  3  3 

stoichiometric  coefficients  v-  and  v.,  the  reverse  reaction  rate 

3  3 

constant  kr>  and  the  equilibrium  constant  K  sk^/k^.  The  rate  constant 
k  is  specified  as  a  function  of  temperature  by  the  expression 
kr  =  AT  exp  (C/T),  and  the  equilibrium  constant  K  by  the  expression 
K  =  A'TB  R(T)  exp  (C'/T).  R(T)  is  the  ratio  of  partition  functions  of  the 
components  on  the  right-hand  side  to  those  on  the  left-hand  side,  given 
in  tabular  form  as  a  function  of  temperature  (see  Category  2  for  table 
parameters).  The  forward  reaction  constant  may  then  be  written: 
k^,  =  AA.'T*B+B  ^R(T)  exp  [{C+C')/T],  Specifications  of  each  reaction 
will  require  the  following  group  of  cards. 


1-40 


41-80 


1-5 


815 


815 


15 


Stoichiometric  coefficients  v.  and 

3 

indexes  of  components  a.,  in  the  order: 


1 


v3  a3 


Stoichiometric  coefficients  v. 


3 


4 

and 


t  * 

indexes  of  components  cu,  as  above 
Temperature  flag  for  K  and  kj..  If 
blank  or  zero,  use  gas  temperature  ; 
if  1  use  electron  temperature 
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Columns 

Format 

6-10 

5x 

May  be  used  for  identification 

11-20 

E10.  3 

Constant  factor  A* 

21-30 

E10.  3 

Temperature  exponent  B* 

31-40 

E10.  3 

Energy  exponent  C' 

41-45 

15 

Temperature  flag  for 

46-50 

5x 

May  be  used  for  identification 

51-60 

E10.  3 

Constant  factor  A 

61-70 

E10.  3 

Temperature  exponent  B 

71-80 

E10.  3 

Energy  exponent  C 

1-80 

8E10.  3 

Table  of  partition  function  ratio 

may  require  more  than  one  card 


Category  4  GEOMETRY  AND  GRID 


1-3 

13 

Number  of  rows 

(max. 

31) 

4-6 

13 

Number  of  columns 

(max. 

81) 

7-9 

13 

Number  of  electrodes 

(max. 

20) 

x-  and  y-spacings  of  the  grid  may  vary  in  x  and  y 
respectively.  In  that  case,  the  number  of  cards  read  will  equal  either 
the  number  of  x-spacings  or  the  number  of  y-spacings,  whichever  is 
greater.  A  uniform  grid  may  be  specified  by  providing  a  single  card 
with  the  negative  of  the  desired  x- spacing  in  the  first  field  and  the 
desired  y-spacing  in  the  second. 


1-10 

11-20 

1.80 


E10.  3 
E10.  3 
2613, 12 


1-80 


2613, 12 


Desired  x-spacing  (m) 

Desired  y-spacing  (m) 

Electrode  geometry  for  the  wall  y  =  0, 
specified  by  giving  the  number  of  grid 
spaces  covered  by  the  leading  insula¬ 
tor,  first  conductor  segment,  next 
insulator  segment,  etc.  More  than 
one  card  may  be  required. 

Electrode  geometry  for  the  wall  y  =  D, 
as  above. 
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Columns  Format 

1-30  3E10.  3  Quadratic  polynomial  coefficients  A,  B,  C 

for  z -dimension  of  channel 
W  =  A  +  Bx  +  Cx2  (m) 

(Needed  for  solution  of  gasdynamic 
problem.  ) 

Category  5  PROBLEM  PARAMETERS 


1-10 

D10.  3 

Total  current  in  x-direction 

(A/m) 

1-80 

8E10.  3 

Array  of  values  of  B,  one  value 

per  grid  column 

(Wb/m2) 

1-80 

8D10.  3 

Total  current  flowing  through 

each  electrode  pair 

(A/m) 

1-80 

8D10.  3 

Given  values  of  streamf unction 

on  each  grid  row  at  column  1 

(if  required) 

(A/m) 

1-80 

8D10.  3 

Given  values  of  streamfunction 

on  last  column  (if  required) 

(A/m) 

1-10 

E10.  3 

■Allowable  relative  variation  in  U 

11-20 

E10.  3 

Allowable  relative  variation  in  T 

21-30 

E10.  3 

Allowable  relative  variation  in  T 

6 

1-80 

8E10.  3 

Allowable  relative  variations  in  < 

:om- 

ponent  number  densities 

1-10 

E10.  3 

Initial  gas  pressure 

(atm) 

11-20 

E10.  3 

Value  of  T  for  initial  current  density 

calculation 

(°K) 

21-30 

E10.  3 

Value  of  T  for  initial  current 
e 

density  calculation 

(°K) 

31-40 

E10.  3 

Seed  fraction  by  weight 

1-10 

E10.  3 

Allowable  relative  error  in  T 

e 

for 

algebraic  solution 

11-13 

13 

Maximum  number  of  iterations  allowed 

in  algebraic  solution  for  Tg 
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Columns 


Format 


Category  6  INITIAL  PROFILES 


1-80 


8E10.  3 


1-80 


1-80 


1-80 


8E10.  3 


8E10.  3 


8E10.  3 


Initial  values  of  all  component 
number  densities  on  rows 
(by  row) 

For  uniform  profile,  make  the 
first  field  of  row  1  specifications 
negative. 

Initial  profile  of  T 


(m“3) 


(°K) 


For  uniform  profile  place  (-T  ) 

© 

in  the  first  field. 

Values  of  T/T^..^  for  all  rows 

For  uniform  profile  place  (-T/T^.^) 
in  the  first  field. 

Values  of  U/U^L  for  all  rows 

For  uniform  profile  place  (-U/U^^) 
in  the  first  field. 


Category  7  AXIAL  VARIATION  OF  U,  T 


1-80 


1-80 


8E1 0.  3 


8E10.  3 


Values  of  T for  all  columns 


(°K) 


For  uniform  distribution  place  (-T^.^) 
in  the  first  field. 

ilues  of  for  all  columns  (m/s) 

For  uniform  distribution  place  (-U^.^) 
in  the  first  field. 


Category  8  ELECTRON  TEMPERATURE  DISTRIBUTION 


1-80 


8E10.  3 


(°K) 


Values  of  electron  temperature  for  all 
grid  points 
First  row-first  column,  first  row- 
second  column,  etc. 
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Columns  Format 

Category  9  RESTART  DECK 

As  suitable  for  each  field  Read  restart  deck  as  punched  by  INLET 


A.  3.  AVAILABLE  OUTPUT 

Both  printed  and  punched  output  are  available  from  INLET, 
the  latter  to  provide  for  continuing  a  calculation  after  interruption  or 
to  facilitate  computer  plotting  of  important  quantities. 


Printed  output,  besides  giving  the  values  at  all  points  of  any 
quantities  described  in  Category  1  of  input,  always  collects  and  prints 
values  of: 

Root -me an- square  variation  of  electron  temperature  between 
present  and  previous  distributions 
Minimum  electron  temperature 
Maximum  electron  temperature 

Number  of  points  where  Tg  is  approaching  a  steady  value  from 
cycle  to  cycle  ("converging*) 

Jy  at  trailing  edge  of  the  middle  electrode  at  y  =  0 
Jy  at  leading  edge  of  the  middle  electrode  at  y  =  D 
Maximum  on  center  line,  and  the  angle  between  the  current 
vector  and  the  y-axis  at  this  point 
Minimum  J^.  on  centerline,  and  the  angle  between  the  current 
vector  and  the  y-axis  at  this  point 
Hall  potential  between  first  and  last  columns 


Each  available  plot  deck  contains  one  card  for  every  grid  point, 
in  the  following  format: 


4-13 

F10.  5 

x- coordinate  of  this  grid  point 
(relative  to  first  column) 

(m) 

14-23 

F10.  5 

y- coordinate  of  this  grid  point 
(relative  to  first  row) 

(m) 

24-33 

F10.  i 

Value  of  chosen  quantity  at  this  grid  point 
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The  restart  deck  contains  values  of  all  the  unknown  quantities 
required  to  continue  a  calculation.  It  should  be  placed  as  a  unit  at  the 
end  of  the  input  deck.  Care  should  be  taken  that  the  input  deck  specifies 
the  same  working  fluid,  geometry  and  grid  as  the  original  computation. 

A.  4.  ROUTINES  USED  IN  PROGRAM  INLET 

Program  INLET  is  composed  of  a  main  program,  which  performs 
control  functions^  and  26  subroutines  which  perform  input,  output  and  the 
tasks  involved  in  the  solution  procedure.  The  list  below  includes  a 
short  description  of  the  function  of  each  routine,  where  M  =  main 
program,  S  =  subroutine,  F  =  function  subroutine. 


ROUTINE 

FUNCTION 

INLET 

M 

Controls  solution 

OP 

S 

Performs  all  printing  output 

IP 

S 

Performs  all  input 

SETUP 

S 

Performs  tasks  associated  with  initializing  or 
continuing  a  problem 

PRPGTE 

S 

Principal  routine  in  solution  of  finite -difference 

equations  for  streamfunction 

SPREAD 

S 

Subordinate  to  PRPGTE 

PSF 

S 

Calculates  values  of  streamfunction  after 

PRPGTE  has  achieved  solution 

CDCF 

S 

Calculates  current  density  field  by  differentiation 

of  streamfunction 

ANODE 

s 

Calculates  elements  of  finite -difference  coef¬ 
ficient  array  along  wall  at  y  =  0 

CORE 

s 

Calculates  elements  of  finite -difference  coef¬ 
ficient  array  away  from  walls 

CATHODE 

s 

Calculates  elements  of  finite -difference  coef¬ 
ficient  array  along  wall  at  y  =  D 

HALL 

s 

Calculates  plasma  properties  for  a  given  set  of 

conditions 

TAND 

s 

Performs  solution  for  gasdynamic  variables 
(U,  T,  Tg,  nQ)  along  streamlines 
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TTEST 

S 

Compiles  and  evaluates  general  characteristics 

POTENT 

S 

Calculates  electric  field  and  potential 

COMPOS 

S 

Calculates  plasma  composition  at  Saha  equili¬ 
brium 

RK 

F 

Calculates  reaction  rate  constants 

DIFTE 

F 

Evaluates  dTg/dx  for  use  in  TAND 

E 

F 

Iterative  solution  of  implicit  algebraic  equation 

for  T 

e 

F 

F 

Iterative  solution  of  implicit  algebraic  equation 

for  T 

e 

ANEWTE 

F 

Subordinate  to  E,  F 

PVAL 

F 

Table  lookup  function 

MATINV 

S 

Matrix  inverter 

PDECK 

S 

Performs  all  punched  output 

PROP 

s 

Calculates  plasma  properties  all  along  a  grid 

row 

RELAX 

s 

Evaluates  relaxation  lengths  and  stabilities 

DZ 

F 

Evaluates  z -dimension  of  channel 

A.  5.  IMPORTANT  VARIABLE  NAMES 

This  list  defines  the  important  variable  names  used  internally 
by  the  FORTRAN  program.  Dimensions  assigned  to  arrays  are  given 
in  parenthesis.  *  =  Double  precision. 

VARIABLE  NAME  DEFINITION 

A(2449,  5)  *  Array  of  finite -difference  coefficients 


AARRAY  (81,  80)  *  Recursion  array  i 

AIXA  *  Total  current  in  x-direction  (A/m) 

AIA(20)  *  Total  current  through  each  electrode  (A/m) 

AJX(31, 81)  Jx  field  (A/m2) 

AJY (31,  81)  Jy  field  (A/m2) 

ALC(20)  Lengths  of  conductor  segments  on 

wall  at  y  =  0  (m) 


-91- 


AEDC-TR-70-86 


VARIABLE  NAME 
ALCC(20) 

ALI(20) 

ALIC(20) 

ALO 

ALOC 

AN(25) 

ANEF(31,  81) 
ANIC(25,  31) 
BARRAY(81,  80)  * 
BF 

BETOSIG(31,  81) 

CP 

CRIT 

D 

EL 

EPSOSIG(3 1,  81) 
EX(31,  81) 

EY  (31,  81) 

HX(81) 

HY(31) 

ICYCLE 

KELECT 

M2 

N 

NCYCLE 

S(2449)  * 


DEFINITION 

Lengths  of  conductor  segments  on 
wall  at  y  =  D 

Lengths  of  insulator  segments  on 
wall  at  y  =  0 

Lengths  of  insulator  segments  on 
wall  at  y  =  D 

Length  of  leading  insulator  on  wall 
at  y  =  0 

Length  of  leading  insulator  on  wall 
at  y  =  D 

Values  of  number  densities  at  the 
current  grid  point 
Electron  number  density  field 
Initial  values  of  number  densities 
Recursion  array  2 
Magnetic  induction 
Ohm's  Law  coefficient  ratio  p/<r 
Specific  heat  c 

p 

Critical  value  of  p/e 
y-dimension  of  channel 
Electron  energy  relaxation  length 
Ohm's  Law  coefficient  ratio  e/<r 
Ex  field 
E  field 

y 

x- spacing  of  grid 
y- spacing  of  grid 
Iteration  cycle  count 
Number  of  electrodes  in  region 
Number  of  grid  rows 
Number  of  grid  columns 
Number  of  iteration  cycles  to  be 
performed 
Streamfunction  field 


(m) 

(m) 

(m) 

(m) 

(m) 

(m-3) 

(m"3) 

(m‘3) 

(Wb/m2) 

(j/(kg  .  °K) ) 

(m) 

(m) 

(v/m) 

(v/m) 

(m) 

(m) 


(A/m) 
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VARIABLE  NAME 
TA(31,8i) 

TE(31,  81) 

U(31,  81) 


DEFINITION 
Gas  temperature  field 
Electron  temperature  field 
Gas  velocity  field 


(°K) 

(°K) 

(m/s) 
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APPENDIX  B 

DESCRIPTION  OF  COMPUTER  PROGRAM  "CHANNEL" 

B.  i.  GENERAL  DESCRIPTION 

Program  CHANNEL  has  been  written  to  solve  the  gasdynamic 
part  of  the  problem  in  the  core  region,  and  to  derive  overall  performance 
data  for  a  J  X  B  device.  Solution  accuracy  is  automatically  held  to 
user's  specifications.  CHANNEL  can  treat  the  flow  of  any  specified 
fluid,  distinguishing  as  many  as  25  separate  components,  formed  from  up 
to  10  distinct  elements  and  taking  part  in  up  to  22  reactions.  Channel 
geometry  and  magnetic  field  and  current  density  distributions  are  required 
input.  Friction  and  heat  losses  to  the  channel  walls,  boundary  layer  and 
electrode  sheath  voltage  drops  may  be  included  at  the  option  of  the  user. 
Printout  of  the  values  of  the  dependent  variables  and  other  important 
parameters  may  be  requested  at  specified  intervals  along  the  channel. 

Upon  reaching  the  specified  terminal  value  of  x,  the  program,  without 
disturbing  the  status  of  the  completed  calculation,  attempts  to  read 
another  input  deck. 

CHANNEL  has  been  coded  in  FORTRAN  for  the  Control  Data 
Corporation  6000  series  computers.  On  a  CDC  6600  computer  the  pro¬ 
gram  occupies  22500  storage  locations,  typically  requiring  40  seconds 
of  central  processor  time  for  a  computation  over  one  meter  of  length  of 
a  J  X  B  accelerator,  with  detailed  output  {at  6  mm  intervals).  Central 
processor  time  is  strongly  affected  by  the  frequency  of  output  required, 
and  by  the  specified  allowable  rates  of  change  of  the  dependent  variables. 

More  detailed  information  relating  to  the  use  of  program  CHANNEL 
is  supplied  in  Sections  B.  2toB.  5  below.  Section  B.  2  contains  a  description 
of  input  data  and  formats,  Section  B.  3  a  list  of  quantities  available  for 
output.  Section  B.  4  contains  a  list  of  the  routines  which  make  up  program 
CHANNEL,  and  Section  B.  5  lists  important  variable  names  used  in  these 
routine  s . 
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B.  2.  INPUT  DATA  AND  FORMATS 
Columns  Format 

1-80  8A10  The  first  card  contains  the  HEADING  to 

be  printed  at  the  top  of  each  new  page. 
Leaving  columns  1-10  blank  terminates 
the  program. 

1-33  1113  The  second  card  is  the  INPUT  CONTROL 

card.  The  remaining  input  is  separated 
into  eleven  categories,  input  of  which 
are  selected  by  punching  the  indexes  of 
those  desired  in  successive  fields  on 
this  card,  beginning  with  the  first  field. 


CATEGORY  INDEX 

Solution  methods  control  1 

Iteration  controls  2 

Specifications  of  geometry, 
current  and  magnetic  field  3 

Fluid  model  4 

Cross  sections  and  weighting 
factors  used  in  Ohm’s  Law  5 

Specification  of  component 


number  densities  to  be  found 


using  differential  equations  6 

Specification  of  chemical 
reactions  7 

Initial  conditions  8 

Electrode  voltage  drops  9 

Required  accuracy  in  integration  10 
Gas  parameters  1 1 


Category  1  METHOD  CONTROL 

1-15  513  Solution  methods  to  be  used  are  selected 

by  punching  their  indexes  in  successive 
fields  of  this  card. 
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Columns  Format 


METHOD  INDEX 

Differential  equation  for  electron 

temperature  1 

Continue  calculation  interrupted 

by  input  2 

Use  as  initial  values  for  all 

variables  those  supplied  in 

the  input  deck  3 

Assume  T  and  all  number 
e 

densities  to  remain  frozen  at 
initial  values  4 

Include  momentum  and  heat 
losses  to  the  channel  walls  5 


Category  2  ITERATION  CONTROL 

1-5  15  Maximum  number  of  iterations  at  any 

point  for  algebraic  solution  of  electron 
temperature  equation. 

6-15  E10.  3  Required  accuracy  in  algebraic  solution 

for  electron  temperature. 

Category  3  GEOMETRY,  CURRENT  AND  MAGNETIC  FIELD 


1-10 

E10.  3 

Initial  value  of  independent  variable,  x 

(m) 

11-20 

E10.  3 

Tabular  interval  (in  x)  for  specification 

of  J  and  J 
x  y 

21-30 

E10.  3 

Terminal  value  of  x 

(m) 

31-40 

E10.  3 

Output  interval 

(m) 

The  y-  and  z -dimensions  of  the  channel 
and  the  magnetic  field  strength  are 
specified  as  quadratic  functions  of  x: 

A  +  Bx  +  Cx^ 
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Columns 

1-30 

31-60 

1-30 


1-80 


1-80 


Format 
3E10.  3 
3E10.  3 
3E10.  3 


8E10.  3 


8E10.  3 


A,  B  and  C  for  y-dimension  (m) 

A,  B  and  C  for  z -dimension  (m) 

A,  B  and  C  for  magnetic  field  strength 

(Wb/m2) 

Values  of  J  and  J  are  specified  for 
x  y 

the  initial  value  of  x  and  at  intervals 
in  x  specified  on  the  first  card  in  this 
category,  up  to  the  terminal  value  of  x. 
Table  of  values  of  J 

y  ? 

(A/rn  ) 


(one  or  more  cards) 
Table  of  values  of  J 

x 

(one  or  more  cards) 


(A/m2) 


Category  4  FLUID  MODEL 


1-5 

6-10 

1-10 

11-20 


15 

15 


Number  of  elements  (up  to  10) 
Number  of  components  (up  to  25) 


21-30 

1-4 

11-60 


Each  element  is  described  on  a  card  of  the  following  format 
A10  Name  of  element 

E10.  3  Fractional  content  of  element  in  the  fluid 

(by  number  density;  required  only  if 
initial  composition  of  fluid  is  not 
specified  in  Category  8) 

E10.  3  Atomic  weight  of  element  (amu) 

Each  component  is  then  described  on  a  card  of  the  following  format 


A4 

1015 


Name  of  component 

The  composition  of  this  component  in 
terms  of  the  specified  elements  is  defined 
by  punching,  in  the  appropriate  fields, 
the  number  of  particles  of  the  correspond¬ 
ing  elements  in  each  particle  of  this 
component 
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Excitation  energy  of  each  component 
above  its  neutral  ground  state.  (eV) 


Category  5  CROSS  SECTIONS 

Cross  sections,  effective  energy-loss  factors  and  weighting 
factors  for  each  neutral  component  are  given  in  tabular  form  as  a  function 
of  temperature.  All  tables  are  specified  by  the  parameters  on  the  first 
card  in  this  category.  The  maximum  table  length  is  51  entries. 

1-10  E10.  3  Minimum  temperature  for  tables  (°K) 

11-20  E10.  3  Tabular  interval  in  temperature  (°K) 

21-30  E10.  3  Maximum  temperature  for  tables  (°K) 

Following  this  card  are  the  tables  of  cross  sections,  energy- 
loss  factors  and  weighting  factors  for  all  neutral  components,  in  the 
order  the  components  were  described  in  Category  4.  The  values  appear 
in  order  of  increasing  temperature, 

1-18  18x  May  be  used  for  identification  or 

sequencing 

19-28  E10.  3  Electron-neutral  collision  cross  section 

Q  (m2) 

29-38  E10.  3  Weighting  factor  A(2) 

39-48  E10.  3  Weighting  factor  A^ 

For  each  neutral  component  there  will  be  up  to  51  cards  in  the 
above  format,  followed  by  the  same  number  in  the  following  format. 

1-10  E10.  3  Effective  energy-loss  factor  for  electron- 

neutral  collisions  5 

eff 


Columns  F  ormat 

1-10  E10.3 
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Columns  Format 

Category  6  INDEPENDENT  COMPONENTS 

1-80  1 61 5  Since  one  may  express  an  algebraic  con¬ 

servation  relation  for  each  element,  it 
is  not  necessary  to  solve  a  differential 
equation  for  every  component.  Those 
components  for  which  differential 
equations  are  to  be  used  are  selected 
by  punching  the  index  of  each  in  succes¬ 
sive  fields,  recalling  the  order  in  which 
they  were  described  in  Category  4. 

This  may  require  two  cards. 


Category  7  CHEMICAL  REACTIONS 


15  Number  of  reactions  (up  to  22) 

Each  reaction,  which  may  be  written  symbolically  as 


Z  via4 

j=l  J 


k'  4 

3*  Z  v!a! 

k  i-i  3  J 
r  3  1 


is  characterized  by  specifying  the  components  a.,  a.  taking  part,  their 

t  J  J 

stoichiometric  coefficients  v.  and  Vj,  the  reverse  reaction  rate 
constant  k  ,  and  the  equilibrium  constant  K=k^/kr>  The  rate  constant 
k  is  specified  as  a  function  of  temperature  by  the  expression 
kr  =  AT  exp  (C/T),  and  the  equilibrium  constant  K  by  the  expression 
K  =  A'T^  R(T)exp  (C'/T).  R(T)  is  the  ratio  of  partition  functions  of  the 
components  on  the  right-hand  side  tc  those  on  the  left-hand  side,  given 
in  tabular  form  as  a  function  of  temperature  (see  Category  5  for  table 
parameters).  The  forward  reaction  constant  may  then  be  written: 
k^  =  AA'T^+^  ^R(T)  exp  [  (C+C')/T] .  Specifications  of  each  reaction 
will  require  the  following  group  of  cards. 
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Columns 

Format 

1-40 

815 

Stoichiometric  coefficients  v-  and 

indexes  of  components  a^,  in  the  order 

V1  ai  V2  a2  V3  a3  V4  fa4 

41-80 

815 

Stoichiometric  coefficients  v.  and 

J 

indexes  of  components  cut,  as  above 
Temperature  flag  for  K  and  k\  If 

1-5 

15 

blank  or  zero,  use  gas  temperature; 
if  1  use  electron  temperature 

6-10 

5x 

May  be  used  for  identification 

11-20 

E10.  3 

Constant  factor  A' 

21-30 

E10.  3 

Temperature  exponent  B* 

31-40 

E10.  3 

Energy  exponent  C' 

41-45 

15 

Temperature  flag  for 

46-50 

5x 

May  be  used  for  identification 

51-60 

E10.  3 

Constant  factor  A 

61-70 

E10.  3 

Temperature  exponent  B 

71-80 

E10.  3 

Energy  exponent  C 

1-80 

8E10.  3 

Table  of  partition  function  ratio  R(T); 

may  require  more  than  one  card 

Category  8  INITIAL  CONDITIONS 

If  Method  3  is  not  selected,  initial  conditions  are  computed  from 
given  stagnation  conditions  and  entrance  Mach  number.  Electron  tem¬ 
perature  is  set  to  a  value  between  Tgta^  and  T  depending  on  a  specified 

parameter  f:  T  =  T  +  f  •  (T  -  T). 

e  siag 

1-10  E10.  3  Stagnation  temperature  Tgta^  {  K) 

11-20  E10.  3  Stagnation  pressure  (atm) 

21-30  E10,  3  Entrance  Mach  number 

31-40  E10.  3  Parameter  f 

If  Method  3  is  selected,  initial  values  of  all  dependent  variables 
must  be  specified. 

1-10  El 0.3  Initial  value  of  velocity  (m/s) 

11-20  E10.  3  Initial  value  of  gas  temperature  (°K) 
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Columns  Format 

21-30  E10.  3  Initial  value  of  electron  temperature  (°K) 

1-80  8E10.  3  Initial  values  of  component  number 

_  3 

densities  (m  } 

May  require  more  than  one  card 

Category  9  ELECTRODE  VOLTAGE  DROP 

The  boundary  layer  voltage  drop  is  specified  as  a  quadratic  poly- 

2 

nomial  in  x,  V  =  A  +  Bx  +  Cx  ,  in  volts 

1-30  3E10. 3  Coefficients  A,  B  and  C  for  boundary 

layer  voltage  drop. 

31-40  E10.  3  Constant  value  of  electrode  sheath  drop  (v) 

Category  10  INTEGRATION  ACCURACY 

Integration  accuracy  is  controlled  by  restricting  the  fractional 
variation  in  every  dependent  variable  over  one  integration  step  to  some 
given  value. 


1-10 

E10.  3 

Allowable  relative  variation  of  velocity 

11-20 

E10.  3 

Allowable  relative  variation  of  temperature 

21-30 

E10.  3 

Allowable  relative  variation  of  electron 

temperature 

31-80 

5E10.  3 

Allowable  relative  variation  of  first  five 

component  number  densities 

1-80 

8E10.  3 

Allowable  relative  variation  in  remaining 
component  number  densities.  May 
require  more  than  one  card. 

Category  11  GAS  PARAMETERS 

1-10 

E10.  3 

Ratio  of  specific  heats  (Cp/cv) 

1-10 

E10.  3 

Recovery  factor 

11-20 

E10.  3 

Wall  temperature  (°K) 

21-30 

E10.  3 

Prandtl  number 
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B.  3.  AVAILABLE  OUTPUT 

In  addition  to  the  current  values  of  all  dependent,,  variables  (u,  T, 
Tg  and  nQ),  printout  at  each  station  includes  the  values  of  the  following 
parameters : 

x-coordinate  of  current  station  Cross-sectional  area  of  channel 
Magnetic  field  strength  Mach  number 

Current  density  components  Mass  density 

Angle  current  vector  makes  with  y-direction 
Gas  pressure  Entropy  relative  to  initial 

Electric  field  components  conditions 

Boundary  layer  voltage  drop  Ohm's  Law  coefficients  cr,  p,  e 

Hall  potential  relative  to  initial  Faraday  potential  E^*  D 
point 

Energy  flux  through  plane  normal  to  flow  direction 
P-MHD  =  f  (E  •  J )  dV  =  Total  MHD  power  input 

P-EL. DROPS  =  f  dx  =  Total  power  dissipated  over  boundary 

layer  voltage  drop 

P- INPUT  s  / (E«7)  dV  +  P-EL.  DROPS  =  Total  power  input  to  flow 
P-WALL  =  /se  dV  =  Total  power  loss  to  walls 
Conversion  efficiency  -  J  (E  •  J  )  dV  /  J  J  Bu  dV 
Total  efficiency  s  f  JyBu  dV/P- INPUT 

B.  4.  ROUTINES  USED  IN  PROGRAM  CHANNEL 

Program  CHANNEL  is  composed  of  a  main  program,  which  per¬ 
forms  control  functions,  and  fourteen  subroutines  which  perform  input, - 
output  and  the  tasks  associated  with  the  numerical  integration.  The  list 
below  includes  a  short  description  of  the  function  of  each  routine, 
where  M  =  main  program,  S  *  subroutine,  F  =  function  subroutine. 

ROUTINE  FUNCTION 

CHANNEL  M  Controls  integration  and  output 

IP  S  Performs  all  input 

OP  S  Performs  all  output 


-102- 


A6DC-TR -70-86 


ROUTINE 

FUNCTION 

SETUP 

S 

Performs  tasks  associated  with  initializing  or 
continuing  a  problem 

STEP 

S 

Performs  Runge-Kutta-Gill  integration  under 
control  of  main  program 

CANON 

S 

Evaluates  components  of  derivative  vector  T 

COMPOS 

S 

Calculates  plasma  composition  at  Saha  equili¬ 
brium  at  a  given  temperature 

RKS 

S 

Calculates  values  of  all  reaction  rate  constants 

at  current  values  of  the  characteristic 

temperatures 

SOURCES 

s 

Evaluates  the  species  source  terms  sq 

SE 

F 

Evaluates  momentum  and  energy  source  terms 

bm  31x3  se 

REMAIN 

S 

Solves  the  system  of  algebraic  equations 
expressing  conservation  of  the  elements 

PROP 

S 

Calculates  electrical  properties  of  the  plasma 

FOX 

s 

Evaluates  parameters  specified  as  functions 

of  x  alone 

EEE 

s 

Performs  iterative  solution  of  implicit  algebraic 
equation  for  electron  temperature 

MATJNV 

s 

Performs  matrix  inversion 

B.  5.  IMPORTANT  VARIABLE  NAMES 


This  list  defines  the  important  variable  names  used  internally  by 
the  FORTRAN  program.  Dimensions  assigned  to  arrays  are  given  in 
parenthesis. 


VARIABLE  NAME 
AMACH 

AN{25)  =  Y  (4-28) 
BETAE 


DEFINITION 

Mach  number 

Component  number  densities 

Ohm*s  law  coefficient  6 

■e 
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VARIABLE  NAME 

DEFINITION 

BETEPS 

Ohm's  law  coefficient  ratio  Pe/ee 

- 

BZ 

Magnetic  field  strength 

(Wb/m2) 

CP 

Specific  heat  c 

P 

(j/(kg.  °K)) 

DELOUT 

Output  interval 

(m) 

ENTHIN 

Enthalpy  flux  at  initial  point 

(J/s) 

ENTROP 

Entropy 

(j/(kg.  °K)) 

EPSE 

Ohm's  law  coefficient  e 

e 

EPSOS 

Ohm's  law  coefficient  ratio  e  / <r 

6 

GAMMA 

Ratio  of  specific  heats  c  /c 

P  v 

HEIGHT 

y-dimension  of  channel 

(m) 

MASFLO 

puA 

(kg/s) 

NCOMP 

Number  of  components 

NELEM 

Number  of  elements 

NREAC 

Number  of  reactions 

A{5)  6 

*  eff’ 

NTAB 

(?\ 

Number  of  table  entries  for  Q,  A'  , 

PIN 

R(T) 

Static  pressure  at  initial  point 

(N/m2) 

PJSQ 

m2 

(A2/m4) 

PJX 

j 

X 

(A/m2) 

PJY 

J 

v 

(A/m2) 

POTION  (25) 

y 

Excitation  energy  of  each  component 

(J) 

POTONS 

Electron  energy  source  ])}e\s. 

(j/(m3-  s) 

PRANDTL 

Prandtl  number 

PRESS 

Static  pressure 

(N/m2) 

PS  TAG 

Stagnation  pressure 

(N/m2) 

RHO 

Mass  density 

(kg/m3) 

SDR  OP 

Sheath  drop 

(v) 

SIGMAE 

Scalar  conductivity 

(mho/m) 

SOURCE  (25) 

Species  source  terms  sq 

T  a  Y (2) 

Gas  temperature 

(°K) 

TE=  Y  (3) 

Electron  temperature 

(°K) 

TOTDEL 

Factor  in  collisional  energy  loss  term:  v.6  ,, 

TSTAG 

Stagnation  temperature 

(°K) 

TW 

Wall  temperature 

(°K) 
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VARIABLE  NAME 

DEFINITION 

U  =  Y(l) 

Gas  velocity 

(m/s) 

VDROP 

Boundary  layer  voltage  drop 

(v) 

WIDTH 

z -dimension  of  channel 

(m) 

XCOOR 

Value  of  x-coordinate 

(m) 

Y  (28) 

Array  of  dependent  variables,  y 

Z(28) 

Array  of  x- derivatives  of  dependent 
variables,  f 
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APPENDIX  C 

DESCRIPTION  OF  COMPUTER  PROGRAM  "LAYER" 

C.  1.  GENERAL  DESCRIPTION 

Program  LAYER  has  been  written  to  solve  the  two-dimensional 
gasdynamic  part  of  the  problem  in  the  compressible,  turbulent  boundary 
layer,  including  magnetohydrodynamic  effects.  LAYER  treats  the  flow 
of  potassium-seeded  air  and  requires  input  only  of  initial  conditions, 
freestream  velocity  distribution,  wall  conditions,  and  electric  current 
and  magnetic  field  distributions.  Upon  reaching  the  specified  terminal 
value  of  x  the  program,  without  disturbing  the  status  of  the  completed 
calculation,  attempts  to  read  another  input  deck. 

LAYER  has  been  coded  in  FORTRAN  for  Control  Data  Corpora¬ 
tion  6000  series  computers.  On  a  CDC  6600  computer  the  program 
occupies  16500  storage  locations.  Because  step  size  in  the  x-direction 
is  limited  to  half  the  boundary i layer  thickness,  typical  execution  time 
requirements  are  difficult  to  state;  however,  each  step  requires  about 
0.  6  second  of  central  processor  time  with  20  grid  points  in  the  y- 
direction,  or  1.  2  seconds  with  40  points. 

More  detailed  information  relating  to  the  use  of  program  LAYER 
is  supplied  in  sections  C.  2  to  C,  5  below.  Section  C.  2  contains  a  descrip¬ 
tion  of  input  data  and  formats,  Section  C.  3  a  list  of  quantities  included 
in  the  output.  Section  C.  4  lists  the  routines  which  make  up  program 
LAYER,  and  Section  C.  5  the  important  variable  names  used  in  these 
routines. 

C.  2.  INPUT  DATA  AND  FORMATS 
Columns  Format 

1-80  8A10  The  first  card  contains  the  HEADING  to 

be  printed  at  the  top  of  each  new  page. 
Leaving  columns  1-10  blank  terminates 
the  program. 


-106- 


AEDC-TR-70-86 


Columns 

Format 

1-5 

15 

If  blank  or  zero,  the  rest  of  this  card 

is  ignored.  Otherwise  the  contents 

of  the  four  remaining  fields  replace 
the  current  values  of  the  correspond¬ 
ing  quantities,  and  no  more  cards 

are  read. 

6-15 

E10.  3 

New  terminal  value  for  x 

(m) 

16-25 

E10.  3 

New  constant  term  in  polynomial 

for  U 

00 

(m/s) 

26-35 

E10.  3 

New  linear  coefficient  in  poly¬ 
nomial  for  U 

oo 

36-45 

E10.  3 

New  quadratic  coefficient  in 

polynomial  for 

1-10 

E10.  3 

Initial  value  of  x 

(m) 

11-20 

E10.  3 

Terminal  value  of  x 

(m) 

21-30 

E10.  3 

Output  interval 

(m) 

31-40 

E10.  3 

Initial  boundary  layer  thickness 

(m) 

41-45 

15 

Desired  number  of  intervals 

across  the  boundary  layer 

1-10 

E10.  3 

Initial  value  of  free -stream 

temperature 

(°K) 

11-20 

E10.  3 

Initial  value  of  static  pressure 

(atm) 

21-30 

E10.  3 

Wall  temperature 

(°K) 

31-40 

E10.  3 

Seed  fraction  (by  weight) 

41-50 

E10.  3 

Effective  energy -loss  factor  in 

electron-neutral  collisions  for 

this  mixture 

1-10 

E10.  3 

Constant  term  in  polynomial 
expression  for 

(m/s) 

11-20 

E10.  3 

Linear  coefficient  in  expression 

for  U 

oo 
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Columns 

21-30 

1-10 

11-20 

21-30 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

1-10 

11-20 

21-30 


Format 
E10.  3 

E10.  3 

EiO.  3 

E10.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 

EIO.  3 


Quadratic  coefficient  in  expression 
for  LT 

oo 

Constant  term  in  polynomial 

expression  for  B  (Wb./m^) 

Linear  coefficient  in  expression 
for  B 

Quadratic  coefficient  in  expres¬ 
sion  for  B 


Constant  term  in  polynomial 

expression  for  Jx  {A/m  ) 

Linear  coefficient  in  expression 

for  J 
x 

Quadratic  coefficient  in  expres¬ 
sion  for  J 

x 

Constant  term  in  polynomial 

expression  for  J  (A/m^) 

Linear  coefficient  in  expression 
for  J 

y 

Quadratic  coefficient  in  expres¬ 
sion  for  J 

y 


Exponent  y  in  expression  for  initial 

profile  of  y:  y^  =  Sfi/N)^  i  =  0.  .  .  N 

for  N+l  grid  points  in  the  y  direction 

Exponent  (3  in  expression  for  initial 

profile  of  gas  velocity: 

U.  =  U  (y./6)P  i  =  0.  ..N 
l  co'Jr 

Exponent  y  in  expression  for  initial 

profile  of  gas  temperature: 

T.  =  T  +  (T  -T  )•  (y./6)v,  i  =  0.  .  .  N 
i  w  oo  w  i 
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Columns  Format 

1-10  E10. 3  Maximum  allowable  relative  error  in 

electron  temperature  calculation 

If  the  exponent  y  for  the  profile  of  y  is  zero,  the  following 
group  of  cards  is  input 

1-80  8E10.  3  Table  of  values  of  y./6  i  =  0.  . .  N 

If  the  exponent  (3  for  the  velocity  profile  is  zero,  the  following 
group  of  cards  is  input 

1-80  8E10.  3  Table  of  values  of  U./U  i  =  0. .  .  N 

1  oo 

If  the  exponent  y  for  the  temperature  profile  is  zero,  the 
following  group  of  cards  is  input 

1-80  8E10.  3  Table  of  values  of  T./T  i  =  0. .  .  N 

1  oo 


C.  3.  AVAILABLE  OUTPUT 


The  values  of  the  normalized  cross -coordinate  go,  as  calculated 
from  the  initial  profiles,  are  printed  on  the  first  page  of  output.  At 
the  initial  station  and  wherever  output  is  requested,  profiles  of  quanti¬ 
ties  dependent  on  y  and  values  of  quantities  independent  of  y  are 


printed. 


Quantities  independent  of  cross -coordinate 


Value  of  x-coordinate  Mass  flow  through  boundary  layer 

Static  pressure,  p  Magnetic  field,  B 

Current  density  components,  Wall  temperature,  T^ 

Jx  and  J  Momentum  thickness,  9 

'  s$c 

Displacement  thickness,  6  x-Reynolds  number,  Re 

Shape  factor,  8  /9 

Displacement  thickness  Reynolds  number,  Reg* 

Momentum  thickness  Reynolds  number,  Reg 

Wall  heat  flux,  JH 

Wall  shear  stress,  t 

'  w 
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1  2 

Skin  friction  coefficient,  c,  =  t  /  r  p  U 

’  f  w  2  roo  oo 

Heat  transfer  coefficient,  cTT  =  JT,  J  p  U  (H  -  H  ) 

*  H  H,  v/  roo  oo  oo  w' 

Quantities  dependent  on  cross -coordinate 


Distance  from  wall,  y 

Gas  temperature  T 

Mass  density,  p 

Mach  number,  M 

Scalar  conductivity,  tr 

Electric  field  components,  E  and  E 

x  y 


Gas  velocity,  u 

Electron  temperature,  Tg 

Electron  concentration,  c 

’  e 

Ionization  relaxation  length, 
Total  enthalpy,  H 


Turbulent  shear  stress,  t 


Electron  number  density,  n 


Ohm's  Law  coefficients,  (3  and  c 
Electron  energy  relaxation  length,  L^ 
Static  enthalpy,  h 


C.  4.  ROUTINES  USED  IN  PROGRAM  LAYER 


Program  LAYER  is  composed  of  a  main  program,  which  con¬ 
trols  the  solution  and  output,  and  subroutines  which  perform  input, 
output  and  various  tasks  associated  with  the  solution.  The  list  below 
includes  a  short  description  of  the  function  of  each  routine,  where 
M  =  main  program,  S  =  subroutine,  and  F  =  function  subroutine. 


ROUTINE 

LAYER 

EEE 

E 

F 

ANEWTE 

BEGIN 

COEFF 

COMPOS 


FUNCTION 

M  Controls  solution  and  output 

S  Controls  iterative  solution  of  implicit  algebraic 
equation  for  Tg 

F  Performs  iterative  solution  for  T 

e 

F  Performs  iterative  solution  for  T 

e 

F  Subordinate  to  E  and  F 

S  Performs  all  input  and  tasks  associated  with 

initiating  or  continuing  computations 

S  Calculates  finite -difference  coefficients  away 

from  the  wall  and  free  stream  boundaries 

S  Calculates  plasma  composition  at  Saha  equilibrium 
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COND 

ENTRN 

LENGTH 

OHM 

OUTPUT 

PAIR 

PAPH 

READY 

SLIP 

SOLVE 
SOLVE  2 

SOURCE 

VEFF 

VISCO 

WALL 

WFi 

WF3 


S  Evaluates  viscosity,  conductivity  and  diffusivity 

factors 

S  Evaluates  mass  entrainment  term 

S  Calculates  boundary  layer  thickness,  ^995*  ant^ 

turbulence  correlation  functions  a^  and  a^ 

S  Calculates  Ohm’s  Law  coefficients  cr,  p,  e,  and 

electric  field  components  E  and  E 

x  y 

S  Performs  all  output 

S  Calculates  static  enthalpy  as  a  function  of 

pressure  and  temperature  for  air 
S  Calculates  temperature  as  a  function  of  pressure 
and  static  enthalpy  for  air 

S  Calculates  distance  from  wall  for  each  grid  point 

S  Calculates  finite -difference  coefficients  for 

points  at  wall  and  free  stream  boundaries 
S  Performs  solution  of  finite -difference  equations 
S  Performs  solution  of  coupled  finite -difference 
equations  for  u,  t 

S  Evaluates  source  terms  for  all  equations 

S  Evaluates  ”eddy  viscosity” 

F  Calculates  molecular  viscosity 

S  Controls  Couette-flow  solutions 

S  Performs  Couette-flow  solution  for  velocity 

S  Performs  Couette-flow  solutions  for  other  unknowns 


C.  5.  IMPORTANT  VARIABLE  NAMES 


This  list  defines  the  important  variable  names  used  internally 


by  the  FORTRAN  program.  Dimensions  assigned  to  arrays  are  given 
in  parenthesis. 


VARIABLE  NAME 

AMACH{43) 

BETAE{43) 

BF 

CF 


DEFINITION 

Mach  number  at  each  grid  point 

Ohm’s  Law  coefficient  P^ 

1  e 

Magnetic  field  strength 
Friction  coefficient  c^ 


(Wb/m2) 


-111- 


AEDC-TR-70-86 


CH 

DD1S 

Heat  transfer  coefficient  c,. 

Displacement  thickness  6 

(m) 

DMOM 

Momentum  thickness  0 

'  (m) 

DX 

Current  interval  in  x  for  calculation 

(m) 

EPSE(43) 

Ohm's  Law  coefficient  e 

e 

EX(43) 

Electric  field  component  E^ 

(v/m) 

EY(43) 

Electric  field  component  E 

(v/m) 

F(4,  43) 

Array  of  unknowns.  F(i,  j)  j  =  1,  43  are 

values  of :  i  =  1  H 

i  =  2  h 

e 

i=3  ce 

i  =  4  T/p 

. 

FLUX  (3) 

Wall  flux  of  H,  h  ,  c 
e  e 

HSTAT(43) 

Static  enthalpy  • 

(J/kg) 

N 

Number  of'intervals  across  boundary  layer 

OM(43) 

Values  of  normalized  coordinate  w 

PEI 

^  Mass  flow  through  boimdary  layer 

(kg/ (m-  s) ) 

.PINF 

Static  pressure 

(.N/m  ; 

RHO{43) 

Mass  density 

(kg/m3) 

SF 

Seed  fraction 

SIGMA(43) 

Scalar  conductivity 

(mho/m) 

TAU(43) 

Turbulent  shear  stress  t 

(kg/ (m*  s  ) 

TEMP(43) 

Gas  temperature  ’  ' 

(°K) 

TEMPE(43) 

1  j  .  "/  .  ’ 

Electron^  temperature 

,(°K) 

TW 

Wall  temperature 

(°K) 

U(43) 

Gas  velocity 

(m/s) 

XNE(43) 

Electron  number  density 

(m-3) 

Y  (43) 

Distance  from  wall 

;  (m) 

YL 

Boundary  layer  thickness  5  <^95 

(m) 
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